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ABSTRACT 


In spite of setbacks due to forest fires, eviction after a change of landowners and unanticipated need to upgrade and 
replace much of the instrumentation, substantial progress has been made during the past three years, resulting in major 
new findings. Although most of the results are in manuscript form, three papers have been published and a fourth was 
recently submitted. 

The data has been subjected to extensive quality control. Extra attention has been devoted to the influence of tilt 
rotation and flux-calculation method, particularly with respect to nocturnal fluxes. 

A relatively complete canopy model (SPA) used in a data assimilation mode can successfully predict the C02 budget 
and NEE by iteratively adjusting the data and model parameters (Chapter 1). To reduce uncertainty in data 
assimilation methods, long-term flux data are needed (both C02 and water vapor exchange). The version of the data 
assimilation method that includes fluxes and state variables is suitable for application in a coupled vegetation- 
atmosphere regional model. 

The SPA model forced by observed meteorological data indicates that in more semi-arid forests, interannual variability 
of precipitation exerts a strong influence on interannual variability of NEE, in contrast to moister regions (Chapter 2). 
These results re-enforce the need to construct a model that correctly couples water vapor exchange and C02 exchange 
in canopy processes and feedbacks between vegetation and the atmosphere. We tested a way to improve the subcanopy 
moisture flux contribution at a boreal site that had a wetter environment (Chapter 3) and found that a representation of 
litter depth in the model could be used to improve subcanopy contributions to energy exchange with the atmosphere. 

Existing canopy models suitable for mesoscale and regional models do not apply to sparse canopies because of the 
importance of daytime buoyancy-generation of mixing in the subcanopy and the controlling influence of subcanopy 
nocturnal inversions, which become much stronger than the stratification above the canopy (Chapter 4). A new 
canopy-mixing model has been developed to accommodate the important influence of subcanopy stability (Chapter 4), 
although further improvements are possible with data taken this year. 

Many Fluxnet sites have reported loss of carbon dioxide inferentially attributed to advection and implying 
accumulation of C02 somewhere downwind. The failure of towers to detect this advective accumulation suggests 
either fundamental observational problems or substantial bias in tower locations. The low-lying young pine site 
appears to accumulate carbon dioxide at night at a rate that exceeds the estimated respiration when mixing and u* are 
weak (Chapter 5), although more investigation is needed. 

Existing techniques for estimating advection of C02 are unreliable, partly due to the inability to estimate the mean 
vertical velocity with sonic anemometers and current flux correction methods (Chapter 6). Estimates of the mean 
vertical motion from the horizontal convergence using a small network of wind sensors appears to be more promising, 
but more expensive in terms of equipment and processing time. More extensive observations will be collected in 
summer 2005. 

Previous/standard methods for calculating nocturnal fluxes with moderate and strong stability are inadequate and lead 
to large random fluxes errors for individual records, due partly to inadvertent inclusion of mesoscale motions that 
strongly contaminant the estimation of fluxes by weak turbulence. Such large errors are serious for process studies 
requiring carbon dioxide fluxes for individual records, but are substantially reduced when averaging fluxes over longer 


periods as in calculation of annual NEE budgets. We have employed a superior method for estimating fluxes in stable 
conditions with a variable averaging width . Mesoscale fluxes are generally unimportant except for events and are 
generally not systematic or predictable. Mesoscale or regional models of our region are not able to reproduce important 
aspects of the diumally varying wind field 


At the writing of this final report, we are working on a short manuscript that re-examines the utility of the commonly 
used u* filter for eliminating records where the NEE is thought to be inadequately estimated due to incomplete mixing. 
We have found that the traditional flux calculation methods for nocturnal turbulence lead to some inaccuracies in the 
u*-filter method. In addition, u* is an incomplete indicator of the ability to assess the degree of mixing since it is 
influenced by pressure fluctuations that do not directly lead to mixing of scalars such as carbon dioxide. Use of the 
heat flux in concert with u* provides a better filter. The net result is that present methods underestimate nocturnal 
respiration. 
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Chapter I 


An improved analysis of forest carbon dynamics using data assimilation 


Global Change Biology (2005) 11, 89-105, doi: 10.1111 /j.l365-2486.2004.00891.x 


An improved analysis of forest carbon dynamics using 
data assimilation 

MATHEW WILLIAMS*, PAUL A. SCHWARZf, BEVERLY E. L AW t, JAMES IRVINEt and 

MEREDITH R. KURPIUSf 

*School of GeoSciences and NERC Centre for Terrestrial Carbon Dynamics , Darwin Building , University of Edinburgh , Edinburgh 

EH9 3JU, UK, f College of Forestry, Oregon State University, Corvallis , OR 97331, USA 

Abstract 

There are two broad approaches to quantifying landscape C dynamics - by measuring 
changes in C stocks over time, or by measuring fluxes of C directly. However, these data 
may be patchy, and have gaps or biases. An alternative approach to generating C budgets 
has been to use process-based models, constructed to simulate the key processes 
involved in C exchange. However, the process of model building is arguably subjective, 
and parameters may be poorly defined. This paper demonstrates why data assimilation 
(DA) techniques - which combine stock and flux observations with a dynamic model - 
improve estimates of, and provide insights into, ecosystem carbon (C) exchanges. We use 
an ensemble Kalman filter (EnKF) to link a series of measurements with a simple box 
model of C transformations. Measurements were collected at a young ponderosa pine 
stand in central Oregon over a 3-year period, and include eddy flux and soil C0 2 efflux 
data, litterfall collections, stem surveys, root and soil cores, and leaf area index data. The 
simple C model is a mass balance model with nine unknown parameters, tracking 
changes in C storage among five pools; foliar, wood and fine root pools in vegetation, and 
also fresh litter and soil organic matter (SOM) plus coarse woody debris pools. We 
nested the EnKF within an optimization routine to generate estimates from the data of 
the unknown parameters and the five initial conditions for the pools. The efficacy of the 
DA process can be judged by comparing the probability distributions of estimates 
produced with the EnKF analysis vs. those produced with reduced data or model 
alone. Using the model alone, estimated net ecosystem exchange of C (NEE) = 
-251±197gCm 2 over the 3 years, compared with an estimate of -419 dr 29 g C m -2 
when all observations were assimilated into the model. The uncertainty on daily 
measurements of NEE via eddy fluxes was estimated at 0.5 g Cm 2 day -1 , but the 
uncertainty on assimilated estimates averaged 0.47 g C m -2 day -1 , and only exceeded 
0.5 g Cm -2 day 1 on days where neither eddy flux nor soil efflux data were available. In 
generating C budgets, the assimilation process reduced the uncertainties associated with 
using data or model alone and the forecasts of NEE were statistically unbiased estimates. 

The results of the analysis emphasize the importance of time series as constraints. 
Occasional, rare measurements of stocks have limited use in constraining the estimates 
of other components of the C cycle. Long time series are particularly crucial for 
improving the analysis of pools with long time constants, such as SOM, woody biomass, 
and woody debris. Long-running forest stem surveys, and tree ring data, offer a rich 
resource that could be assimilated to provide an important constraint on C cycling of 
slow pools. For extending estimates of NEE across regions, DA can play a further 
important role, by assimilating remote-sensing data into the analysis of C cycles. 

We show, via sensitivity analysis, how assimilating an estimate of photosynthesis - 
which might be provided indirectly by remotely sensed data - improves the analysis 
of NEE. 


Correspondence: Mathew Williams, fax 4- 44 131 662 0478, 
e-mail: mat.williams@ed.ac.uk 
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estimates of ecosystem C stocks and fluxes, with 
reduced uncertainty compared with the original ob- 
servations, or the model alone. The argument of this 
paper is that combining measurements and modelling 
through DA generates more precise estimates of C 
dynamics, and simultaneously highlights areas where 
model improvement is required. 

Methods 

The premise of DA is that neither models nor 
observations can perfectly describe a system, but an 
analysis that combines model and data will provide a 
better estimate of system dynamics than model or 
observations alone. DA is a process for the optimal 
combination of information about a system, which 
evolved from the engineering approaches to filtering 
and control theory applied in missile guidance and 
interception (Maybeck, 1979). DA has been applied in 
meteorology for forecasting (Lorenc, 1986), and ex- 
tended into stream ecology (Cosby, 1984), oceanogra- 
phy (Eknes & Evensen, 2002), and soil science 
(Heuvelink & Webster, 2001) for time series analysis. 

DA is the process of finding the model representation 
that is most consistent with observations (Lorenc, 1995). 
DA recognizes that there are never sufficient observa- 
tions to represent the state of a system at any one time. 
For a detailed, complete picture, further information is 
required, such as knowledge of the behaviour and 
probable structure of the system. In DA, knowledge of 
system evolution in time is usually embodied in a 
model. In sequential assimilation, the approach we 
demonstrate here, the model organizes and propagates 
forward information from previous observations 
(Lorenc, 1995). When new information becomes avail- 
able, the prediction, or forecast, of the model can be 
compared with these observations and corrected. A poor 
model will drift and will be frequently and heavily 
corrected; an effective model will require little reinitializa- 
tion by observations. However, it is not simply a question 
of fitting the model to the new data, as the assimilation 
process must also conserve the information provided by 
the model itself and by previous observations. 

The DA technique that we use here is the Kalman 
filter (KF) (Kalman, 1960), which has been widely used 
(Grewal, 1993), and, given various assumptions, has 
been shown to be an optimal, variance-minimizing 
analysis (Maybeck, 1979). The basic KF requires three 
assumptions: that a linear model can describe the 
system, and that the system and measurement noise are 
both white and Gaussian. Developments of the basic KF 
have provided means to deal with deficiencies in these 
assumptions (Grewal, 1993). The product of the KF is 
an estimate, or analysis, of the state variables that takes 


account of prior knowledge plus new observations to 
ensure that estimated errors are statistically minimized. 

The KF generates a variance-minimizing analysis by 
combining the model forecast with the observations, 
weighted according to these prediction and measure- 
ment error covariances. If measurement noise is large, 
relatively less emphasis is placed on the current 
observation. If measurement noise is small, or the 
model estimation error is large, the corrected estimate 
approaches the observation. As inputs, the KF requires 
details on the covariances of both the model forecast 
and of the measurements. The error covariance of the 
measurements can be derived from knowledge of the 
accuracy of the techniques used to measure them. The 
error covariances for the model prediction are com- 
puted by solving an equation for the evolution in time 
of the error covariance matrix of the model state 
(Evensen, 2003). 

The KF uses covariance matrices to store information 
on the uncertainty in the models, observations, and 
analysis. Storing, integrating, and inverting large 
covariance matrices is computationally expensive and 
the matrix operations required for the analysis are not 
always robust. Evensen (1994) suggested that, instead 
of storing a full covariance matrix, the same error 
statistics can be represented approximately using an 
appropriate ensemble of model states. The ensemble KF 
(EnKF, Evensen, 1994) uses a Markov Chain Monte 
Carlo method to solve the time evolution of the 
probability density of the model state. The probability 
density is represented by a large (ca. 100-1000) 
ensemble of model states, and these are integrated 
forward in time by a differential equation (i.e., the 
model) with a stochastic forcing term representing the 
model errors. Each ensemble member evolves in time 
according to 

^fy +1 = M(^y) + dcjj , (1) 

where ip is the state vector, j counts from 1 to N, where 
N denotes the number of model state ensemble 
members, k denotes time step, M is the model operator 
or transition matrix, and dq is the stochastic forcing 
representing model errors from a distribution with 
mean zero and covariance Q (Burgers et al ., 1998). 

Similarly, observations are treated as random vari- 
ables by generating an ensemble of observations from a 
distribution with the mean equal to the measured value 
and a covariance equal to the estimated measurement 
error (Burgers et al., 1998). Thus, we define the new 
observations 

dj = d + Ej , (2) 

where d are the observations, and e are drawn from a 
distribution of zero mean and covariance equal to the 
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Introduction 

Quantifying the carbon (C) dynamics of the terrestrial 
biosphere is a current and major concern for earth 
system science. A major uncertainty concerns identify- 
ing whether regions or landscapes are sources or sinks 
for C. There are two broad approaches to quantifying 
landscape C dynamics - by measuring changes in C 
stocks over time, or by measuring fluxes of C directly. 
Stock inventories involve recording the mass of C in 
living biomass (leaves, stems, and roots), and in litter 
and soil pools, and quantifying the changes in pool 
sizes between sampling times (Turner et al, 1995; Malhi 
et al, 2002). The advantage of this approach is that 
measurements are generally cheap and simple, 
although labour intensive. The difficulty is that pools 
are spatially patchy (e.g., a large proportion of C can be 
in a few large tree stems), belowground pools (such as 
fine roots, root litter, soil C) are difficult to measure, and 
that monitoring is generally restricted to small plots by 
logistical constraints. 

Direct determination of C fluxes has been revolutio- 
nized by the development of automated measurement 
systems of net carbon dioxide (C0 2 ) and water vapour 
exchange between land and atmosphere, via the eddy 
covariance technique (Baldocchi, 2003). Cuvettes and 
chamber measurements of soil effluxes, stem and 
foliage respiration, and leaf photosynthesis have 
enhanced the understanding of processes contributing 
to the eddy fluxes (Law et al, 1999a, 2001a). Flux data 
can be generated somewhat continuously (e.g., half- 
hourly), and the instruments sample a 'footprint' of the 
surrounding landscape covering a few km 2 . However, 
these data often have gaps, and filling these may 
introduce biases, and increase uncertainty. Also, night- 
time flux data can be biased when winds are light and 
intermittent (Goulden et al, 1996), and complex terrain 
may jeopardize some of the assumptions of the 
approach (Finnigan et al, 2003). Finally, the scale of 
measurement is often uncertain because the footprint 
varies depending on wind speed and direction (Schmid 
& Lloyd, 1999). 

An alternative approach to generating C budgets has 
been to use process-based models, constructed to 
simulate the key processes involved in C exchange 
(Farquhar & von Caemmerer, 1982; Jarvis et al, 1985; 
Parton et al, 1988). The advantage of using models is 
that they can be extended across large spatial domains 


and into the future, given the relevant driving variables 
(Running et al, 1999; White et al, 2000; Rastetter et al, 
2003). Forecasts are possible because models incorpo- 
rate a representation of the simulated system and its 
dynamics. Rules, such as the conservation of mass, can 
be enforced to guide system trajectories. The disadvan- 
tage of modelling is that the process of model 
construction is arguably subjective. Occam's razor - 
making models as simple as possible, but no simpler - 
is a useful guiding principle. But there is always a 
danger that the model's representation of the system is 
not accurate. Other problems include parametrization - 
setting the rate constants on fluxes in a compartment- 
flow model, for instance. Generally these parameters 
are unknown and have to be derived from data, 
somehow. There is always a danger that poorly defined 
parameters will be 'tuned' to give good output. When 
several parameters are tuned, the right answer may be 
generated for the wrong reason (Williams et al, 2001a). 

Generally, scientific papers on C budgets can be 
classified into stock (Phillips et al, 1998), flux (Wofsy 
et al, 1993), or model approaches (McGuire et al, 1993), 
following the definitions given above. However, some 
papers do attempt to combine the methods. Ecological 
inventories are now being compared with microme- 
teorological data (Ehman et al, 2002), and models are 
often corroborated against observations (Running, 1994; 
Williams et al, 1996; Law et al, 2003). 

Models are generally parametrized with some subset 
of observational data, and tested against remaining 
data. Such tests are designed to show that the model 
can effectively describe the observed system by 
demonstrating a strong correlation, or a low mean 
error, between prediction and observation. This stan- 
dard modelling approach assumes a primacy of the 
data over the model representation. But if this is the 
case, then the standard approach is inefficient. It would 
be much more worthwhile to use all the available data 
to improve the model and minimize confidence limits 
on predictions, rather than only a subset. Here we 
demonstrate the technique of data assimilation (DA), 
combining all available data with a model, to develop a 
full analysis of C cycling in a ponderosa pine forest. 

Objectives 

The objective of this paper is to demonstrate why the 
application of DA techniques results in improved 
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estimated measurement error. The analysis step of the 
EnKF updates each of the model state ensemble 
members using the following equation: 

= + (3) 

where ip* is the forecast state vector (i.e., the prediction) 
and ip a is the analysed estimate generated by the 
correction of the forecast. H is the observation operator, 
a matrix that relates the model state vector to the data, 
so that the true model state is related to the true 
observations by 

d ' = Hip*. (4) 

K<> is the KF gain matrix, which determines the 
weighting applied to the correction (Burgers et al, 
1998). 

We use the EnKF based on its ability to predict error 
statistics for strongly nonlinear systems and for its 
simplicity and numerical stability (Evensen, 2003). The 
analysis generated by the EnKF is a combination of the 
model forecast and a number of influence functions, 
one for each of the measurements (Allen et al, 2002). 
These influence functions are computed from the 
ensemble statistics, and thus include cross-correlations 
between the different variables in the model. The 
influence functions summarize the correlations be- 
tween model state variables that are determined from 
the ensemble of model runs. In effect, the ensemble 
serves as a large sensitivity analysis, quantifying how 
small changes in each state variable affect system 
dynamics. Changes to one variable resulting from a 
new observation are transmitted to other model 
variables according to the strength of these cross- 
correlations. This is a particularly powerful component 
of the EnKF - information on a single observation is 
transmitted to all state variables via the connections set 
out in the model. The FORTRAN code required to 
perform the EnKF is provided in Evensen (2003). 

The study area 

The Metolius young ponderosa pine site is located in a 
Research Natural Area (44°26'N, 121°34 / W, elevation 
1165 m) in the eastern Cascades, near Sisters, Oregon. 
The site was clear cut in 1978, and has regenerated 
naturally since then, with some recent thinning. In 2002, 
there were 431 trees ha~\ with a mean diameter at 
breast height (DBH) of 11.3 cm. The understorey 
vegetation is sparse with patches of bitterbrush ( Purshia 
tridentata) and bracken fern ( Pteridium aquilinum), and a 
groundcover of strawberry ( Fragaria vesca). The site is in 
a semiarid region that experiences warm, dry summers 
and wet, cool winters. While the flux tower samples a 
variable footprint, most of the remaining data were 


collected within a lOOmxlOOm plot located ~ 50m 
upwind of the tower. 

Flux measurements 

A variety of flux data were used to generate daily 
estimates of specific C fluxes for use in the analysis. 
We made continuous eddy covariance measurements 
to determine half-hourly fluxes of C0 2 throughout 
2000-2002 at 12 m height, ~ 9 m above the canopy 
(Law et al, 1999a, b; Anthoni et al, 2002). Data were 
screened to remove possible eddy covariance instru- 
mentation and sampling problems (e.g., friction velo- 
city, u*<0.2ms" 1 ), and fluxes were also rejected when 
unreasonably large C0 2 fluxes ( I F c I >25 |imol m -2 s _1 ) 
were observed. C0 2 concentration was measured every 
30 min at 1, 3, and 12 m above the ground. The buildup 
or release of C0 2 from within the canopy was 
quantified by determining the rate of change of C0 2 
along the vertical profile. The total C0 2 flux was then 
calculated as the sum of C0 2 flux and this change in 
storage. Data gaps were filled based on seasonal 
empirical relationship with environmental variables 
(PAR, VPD) derived from valid flux data (Anthoni et al, 
1999). We generated daily net ecosystem exchange of 
C0 2 (NEE) data for days in which less than 25% of the 
48 possible half-hour measurements were gap-filled; for 
the 3-year period of this study, this amounted to 684 
daily NEE values. 

Other periodic gas exchange measurements included 
foliage and stem respiration, and soil surface C0 2 
fluxes (Law et al, 1999b). Soil respiration was measured 
using six automated chambers installed in 1999 (Irvine 
& Law, 2002); total daily effluxes were recorded on 401 
days during 2000-2002. Root contributions to soil C 
effluxes were measured by recording C0 2 efflux with a 
manual system, removing a 30 cm soil core, extracting 
the roots, and then measuring root C0 2 effluxes 
(methods in Law et al, 2001a). These soil/root data 
were collected on 3 days in 2000 (days 153, 201, and 
278), and 2 days in 2001 (115 and 227). Pine foliage 
respiration was determined using a temperature func- 
tion fitted to 12 nights of cuvette data collected during 
1999 and 2001 (Law et al, 1999b); shrub foliage 
respiration rates were determined similarly using five 
nights of data (Law et al, 2001a). Site-level estimates of 
total foliar respiration were determined from air 
temperature data, temperature response functions for 
trees and shrubs, and partitioning of estimated site leaf 
area index (LAI) between trees and shrubs. We 
estimated sapwood respiration at the Y site using a 
temperature response function developed from data 
at a nearby old-growth site, using sapwood volume 
estimates derived from tree rings at the young site 
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(Law et al, 1999a, 2001a). Because sapwood of young 
trees has higher respiration rates than old trees, we 
likely underestimated the respiration loss from the 
system, but in the old forest, sapwood respiration 
accounted for only 10% of total ecosystem respiration 
(Law et al . , 1999a). 

We limited estimates of total ecosystem respiration to 
those 401 days where soil measurements were avail- 
able. Estimating foliar and stem respiration on these 
days required temperature and biomass data for each of 
these days. Errors will be dominated by uncertainty in 
the foliar biomass and sapwood volume estimates, 
although the latter should be relatively small. To 
estimate total autotrophic respiration, we combined 
estimates of the root fraction of soil effluxes with foliar 
and stem respiration measurements. Because there 
were no automated soil chamber data within 55 days 
of the root respiration measurement on day 115 in 2001, 
we discarded this data point. Thus, over the 3-year 
period there were 4 days when we were able to produce 
estimates of autotrophic respiration. 

To provide independent estimates of gross primary 
production (GPP) (direct observations were lacking), we 
used the soil-plant-atmosphere (SPA) model (Williams 
et al, 1996), which is based on the underlying biochem- 
istry of carboxylation, of light interception, gaseous C0 2 
exchange, and the impacts of soil moisture on stomatal 
conductance. The SPA model has multiple canopy layers 
and a 30 min time step, and has been previously 
parametrized and applied in ponderosa pine ecosystems 
(Williams et al, 2001b), generating reasonable estimates 
of C fixation throughout the annual cycle. The SPA model 
requires as driving variables a daily phenology of LAI, 
foliar N, and root biomass, which we generated from 
interpolations of the limited data available over time. We 
used measured estimates of maximum carboxylation and 
electron transport rates on pine foliage at the site to 
parametrize the model (Law et al . , 2003), and time series 
of sap flow data to parametrize the hydraulic character- 
istics of the trees (Schwarz et al, in press). The close 
correspondence of SPA predictions of daily transpiration 
with sap flow estimates ( R 2 varied from 0.82 to 0.87 over 
the three years) means that we have confidence in the 
SPA predictions of stomatal opening (Schwarz et al, 
2004). The SPA model, in effect, translates sap flow 
observations into GPP estimates, which are then assimi- 
lated. This means that NEE data, respiration data, and 
GPP pseudo-observations can all be assimilated into the 
analysis and their consistency can be determined. 

Stock measurements 

We estimated LAI (one-sided LAI) from optical 
measurements with an LAI-2000 plant canopy analyzer 


(LI-COR, Lincoln, NE, USA), on a 10 m grid within the 
100 m x 100 m plot. We also measured needle clumping 
within shoot, and clumping at scales larger than shoot 
to correct for these effects on LAI estimates. LAI data 
were collected at four times through the 3-year period 
of this study - once in 2000 and in 2001, and twice in 
2002. Foliar mass was determined from LAI and direct 
measurements of C mass per unit leaf area on foliage 
samples. 

On five 10 m radius subplots within the main plot, 
we recorded dimensions of all trees (stem height, and 
DBH, 1.37 m) and shrubs (length, width, height, 
and diameter at shrub base). Aboveground biomass 
and coarse root mass were then determined from 
allometric relationships derived from destructive har- 
vest of five trees covering the range of sizes present 
(Law et al, 2001b). Likewise, shrub biomass was 
determined from the destructive harvest of five to nine 
shrubs per species and scaled to the site by the number 
of shrubs per size class (Law et al, 2001b). Dimensions 
were sampled at three times during the 3-year period of 
the study, once in 2000 and twice in 2002. 

Litterfall production (<lcm diameter) was deter- 
mined from monthly collections of litterfall in 20 trays 
(0.13 m 2 each); litter was separated into foliage and 
woody material (Law et al, 2001b). Litterfall was only 
collected during the latter half of each of the 3 years, 
when the majority of litterfall occurs, resulting in 
around 18 months of data over 3 years. We estimated 
daily litterfall for both foliage and woody material on 
the basis of a constant rate across the month. 

Fine root (<2mm diameter) biomass was estimated 
from soil cores collected on 2 days during 2002 (days 
134 and 219). On each occasion, cores were extracted 
from 18 different locations (six sampling points along 
three transects) using a 7.2 cm diameter corer. At each 
sampling location, samples were divided into three 
layers: 0-20, 20-50, and 50-100 cm depth. All samples 
were refrigerated until processed in the laboratory, at 
which time the roots were separated from the 
surrounding soil and sorted according to diameter 
class and then further separated into live and dead 
categories based on visual inspection of each root 
segment. Following the separation and sorting, the root 
samples were dried at 70 °C for 48 h and then weighed 
to measure biomass. 

C storage in soil was determined from 27 soil cores 
to 1 m depth. Live vegetation and roots were removed 
and total C was determined to 1 m depth (Law et al, 
2003). Coarse woody debris (CWD) was measured 
in four 75 m line transects, and fine wood debris 
was recorded on four 15 m transects per subplot for 
all four subplots, following the protocol of Harmon & 
Sexton (1996). Six samples of forest floor fine litters 
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Table 1 Estimated errors associated with each component of the model and with each set of observations, and the total number of 
daily observations available 


Symbol 

Model error 

Observational 

error 

Number of 
observations 

Pool/ flux description 

c, 

20% 

10% 

4 

Foliar C mass 

C„ 

20% 

10% 

3 

Wood C mass 

C r 

20% 

30% 

2 

Fine root C mass 

C„, 

20% 

30% 

1 

Fresh litter C mass 

CsOM/WD 

20% 

30% 

1 

Soil organic matter plus woody 
debris C mass 

A, 

20% 

N/A 


Foliage allocation rate 

Ay, 

20% 

N/A 


Wood allocation rate 

a t 

20% 

N/A 


Fine root allocation rate 

U 

0.5 g Cm -2 day -1 

20% 

18 

Foliage litter production rate 

Lw 

0.5 g Cm -2 day -1 

20% 

18 

Wood litter production rate 

Lr 

0.5 g Cm -2 day -1 

N/A 


Fine root litter production rate 

D 

20% 

N/A 


Decomposition rate of litter 

G 

20% 

30% 

1096 

Gross primary production 

Ktot 


20% 

401 

Total respiration rate 


20% 

50% 

4 

Autotrophic respiration rate 

^hl 

20% 

N/A 


Heterotrophic respiration rate of fresh litter 

&h2 

20% 

N/A 


Heterotrophic respiration rate of SOM/WD 

NEE 

N/A 

0.5 g Cm -2 day -1 

684 

Net ecosystem exchange rate of C 


For each time step, the predictive error on each pool and flux is drawn from a normal distribution of mean zero and with a standard 
deviation defined as a percentage of the mean pool or flux at that time. For modelled litter fluxes and observed net ecosystem 
exchange (NEE), an absolute and constant value is set instead (see text). 


were collected to determine nonwoody litter mass (Law 
et al . , 2003). 

Confidence limits on observations 

Correct estimation of observational error is crucial to 
the quality of the analysis (Table 1), because the 
magnitude of observational errors determines to what 
extent the simulated fields will be corrected to match 
the observations. Observation error variances are best 
specified by knowledge of instrumental characteristics, 
and by comparing replicated samples. Observation 
error correlations are usually assumed to be zero - 
distinct measurements are assumed to be affected by 
physically independent errors. In most cases, we 
specify the standard deviation of model error as a 
percentage of the mean. 

Errors in the total respiration estimates arise from the 
limited sampling of the automated soil chambers in 
space, and the use of empirical relationships to scale 
foliage and stem respiration in time. Because soil 
respiration is the largest component of ecosystem 
respiration, we set the coefficient of variation in total 
respiration in accordance with the reported values for 
the soil chamber data at the study site, ~ 0.2 (Irvine & 
Law, 2002). The uncertainty of autotrophic respiration 


estimates must be relatively larger than that of total 
respiration, because the autotrophy is estimated by a 
highly invasive approach, so errors are set to 50%. A 
detailed uncertainty analysis of modelling of GPP 
suggests errors of up to 30% (Williams et al. , 2001c), 
and this value was assigned to SPA model estimates of 
GPP. An analysis of the uncertainties in calculating 
daily NEE from eddy covariance (Anthoni et al., 1999) 
suggests an error on daily estimates of 12%. The 
uncertainty because of potential advection at night 
and during early morning and late evening transition 
periods could not be assessed. In addition, recently 
discovered software error in the LI7500 flux instrument 
suggests that the flux estimates of NEE are biased high 
by about 20% (biased towards a stronger sink for C0 2 ; 
LI-COR Inc.). However, the nature of NEE (it can be 
positive or negative) means that defining errors by a 
coefficient of variation is unsuitable, so, instead, errors 
are set at 0.5 g Cm -2 day -1 . 

Errors on pool measurements (Table 1) were deter- 
mined from the variation in replicated measurements. 
The coefficient of variation on LAI observations at the 
site was 10%, and stem sampling indicated a standard 
error on wood biomass equal to 10% of the mean (Law 
et al., 2001b). Based on replicated soil cores, we estimate 
the error on fine root biomass at 30% (data not shown). 
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Models 

Description of C dynamic model 

We represent the C cycle with a simple box model of 
pools connected via fluxes (Fig. 1). There are five pools 
- C content of foliage (Cf), woody stems and coarse 
roots (C w ) and fine roots (C r ), and of fresh leaf and fine 
root litter (Cutter) and soil organic matter (SOM) plus 
WD (Csom/wd)- There is one pseudopool, the daily 
accumulation of photosynthate (GPP) that is entirely 
utilized each day. The fluxes among pools are based on 
the following assumptions: 

1. All C fixed during a day is either expended in 
autotrophic respiration or else allocated to one of 
three plant tissue pools - foliage, wood, or fine roots. 

2. Autotrophic respiration is a fixed fraction of total 
photosynthetic fixation (Thornley & Cannell, 2000), 
and it is not directly temperature sensitive. 

3. Plant allocation and litterfall are donor-controlled 
functions with no direct environmental influence 
and constant rate parameters. 

4. Soil transformations are sensitive to temperature, 
with a Q 10 of 2.0. Otherwise, the only environmental 
forcing in the C model is on GPP, via solar radiation, 
air temperature, and soil moisture. 

5. All C losses are via mineralization; there is no 
dissolved loss term. 

The aggregated canopy model (ACM) of photosynth- 
esis (Williams et al, 1997) provides the forecast estimate 
of daily C inputs to the system. The ACM is a big-leaf, 
daily time step model that estimates GPP as a function 
of LAI, foliar nitrogen, total daily irradiance, maximum 
and minimum daily temperature, day length, atmo- 
spheric C0 2 concentration, soil-plant water potential, 
and total soil-plant hydraulic resistance. The model has 
10 parameters, and we use a local model calibration 
(see Appendix). Using measurements of leaf C mass per 


area (111 g Cm -2 leaf area), estimates of LAI for the 
ACM are determined directly from the C { pool. For 
simplicity we assume a constant value for foliar N 
concentration (2.7 g N m -2 leaf area). To characterize the 
changes in moisture limitation through the summer, we 
drive the ACM model with estimates of soil-leaf water 
potential difference and hydraulic resistance, both 
derived from a detailed ecophysiological study under- 
taken at the site using the SPA model (Schwarz et al, 
2004). 


Initial conditions and rate constants 

The C box model has nine unknown parameter 
constants (Table 2), and the initial values of the five C 
pools are also unknown (Table 3). To generate estimates 
of parameters and initial conditions, we nested the EnKF 
within an optimization routine. This routine varied the 
unknown parameters and initial conditions to find the 



Fig. 1 C dynamic model, showing pools (boxes) and fluxes 
(arrows). Allocation fluxes are marked A, litterfall fluxes by L, 
and respiration by R (split between autotrophic (a) and 
heterotrophic (h)). D is a decomposition flux and GPP is gross 
primary production. Temperature-controlled fluxes are marked 
by dashed lines. 


Table 2 Estimated values of the model parameters as calculated by the optimization routine 


Parameter 

Description 

Model value 

Low /high bounds 

h 

Decomposition rate constant 

4.4 x 1(T 6 

1 X 10' 6 /0.01 

*2 

Autotrophic respiration as a fraction of GPP 

0.47 

0.2/0.7 

*3 

Fraction of NPP allocated to foliage 

0.31 

0.01/0.5 

U 

Fraction of NPP allocated to fine roots 

0.43 

0.01/0.5 

t 5 

Turnover rate of foliage 

2.7 x 10“ 3 

2 x 10~ 4 /0.02 

*6 

Turnover rate of woody matter 

2.06 x 10“ 6 

2 x 10~ 6 /0.02 

t7 

Turnover rate of fine roots 

2.48 x 10' 3 

2 x 10“ 4 /0.02 

t8 

Mineralization rate of fresh litter 

2.28 x 10' 2 

5 x 10 -s /0.5 

*9 

Mineralization rate of soil organic matter and woody debris 

2.65 x 10' 6 

1 x 10 _6 /0.5 


Low and high bounds for each parameter, which define the search space are also shown. 
NPP, net primary production; GPP, gross primary production. 
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Table 3 Estimated values of the initial conditions in each dynamic C pool 


Symbol 

C pool 

Initial value (g C m 2 ) 

Low /high bounds 

C, 

Foliage 

58 

30/70 

C w 

Wood (stems and coarse roots) 

770 

750/900 

C r 

Fine roots 

102 

50/300 

C lit 

Fresh foliar and fine root litter 

40 

20/3000 

CsOM/WD 

Soil organic matter plus woody debris 

9897 

6000/10000 


Low and high bounds for each pool, which define the search space. 


values that minimized the sum-of-squared differences of 
the innovations (i.e., the difference between model 
forecast /prediction and observations) for all available 
observations (Table 1). This approach, in effect, under- 
takes numerous implementations of the EnKF with 
varied initial values and parameters, and identifies in 
which implementation the predictions require the mini- 
mum correction. This implementation is then assumed 
to have the optimal parameter set and initial conditions. 

The optimization routine uses the quasi-Newton 
method and a finite difference gradient, UMINF, IMSL 
Math Library. We set specific bounds to the parameters 
and initial conditions to ensure realistic values (Tables 2 
and 3). We used the optimized parameters in all the 
subsequent analyses, although we also undertook a 
sensitivity analysis to find how variation in parameters 
affected the analysis, and also affect model forecast bias. 

Model error variances and ensemble size 

The model errors determine how rapidly the uncer- 
tainty on the forecast of the state vector grows over 
time. The errors are determined by random distur- 
bances to the system, errors in measured inputs 
(drivers), and the error inherent in representing a 
complex system as a simple model (Cosby et al ., 1984). 
The standard deviations of model error can be 
expressed as a fraction of mean values, and this is 
useful as it largely restricts the ensembles to positive, 
and thus meaningful, values when the mean of the 
particular ensemble is close to zero (e.g., GPP during 
winter). We set the model uncertainties to 20% (Table 1), 
and undertook sensitivity analyses to explore the 
implications of changes in these values. As an excep- 
tion, we assigned a relatively large, constant error to 
litterfall simulations (Table 1). The modelling of litter- 
fall uses a simple rate constant, and so does not 
simulate stochastic events, such as storms, or periodic 
events, such as leaf senescence, that cause pulses of 
litter. By setting a high and constant error on litterfall 
modelling, any available observations are used to form 
the majority of the analysis, and we are explicitly 
recognizing the weakness of our litterfall submodel. 


The ensemble size of the EnKF is the number of 
model states that are concurrently predicted and 
analysed. The differences between the analyses form a 
statistical sample of analysis errors. We set ensemble 
size to 200, large enough to ensure the correct estimate 
of the error variances in the predicted model state 
(Allen et al , 2002). 


Results 

Parameter estimates and initial conditions 

The fitting process generated estimates of the nine 
model parameters and five initial conditions using all 
the available flux and stock data. The optimization 
estimated that 47% of GPP was respired by plants each 
day (Table 2), and 53% remained for net primary 
production (NPP). Of NPP, 31% was allocated to foliage 
each day, 43% to fine roots, and 25% to woody stems 
and coarse roots. The parametrized daily turnover rates 
of plant tissues (Table 2) suggest a leaf life span of 1 .0 
years, a fine root life span of 1.1 years, and a life span of 
woody materials of 1323 years. The parametrized 
estimate of fresh fine litter turnover, via mineralization 
and decomposition combined, was 0.1 years at a 
constant 10 °C. Parametrized turnover of SOM and 
WD was 1033 years, at a constant 10 °C. 

The estimated initial conditions (Table 3) assigned a 
smaller C mass to foliage than to fine roots, but 83% of 
the initial vegetation C pool was allocated to woody 
biomass. The C mass of CWD and SOM are more than 
an order of magnitude greater than total plant C 
biomass. 


Analysis of carbon ecosystem exchanges 

Over the 3 years, 2000-2002, the analysis suggests a 
total NEE of -419 ± 29 g Cm" 2 (i.e., a clear C sink) 
(Table 4). The total estimated GPP over the 3 years was 
2172 ± 18 g Cm -2 , with 2002 being the most productive 
year (Table 4). The analysis suggests that autotrophic 
sources contributed to 58% of total respiration (Table 4) 
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Table 4 Analysis of annual C fluxes (gCm 2 yr ? ) for 3 years, and their sum total 


Flux 

2000 

2001 

2002 

Total 

SD 

R 2 

N 

RMSE 

Gross primary production 

716 

702 

753 

2172 

18 

0.95 

1096 

0.32 

Total respiration 

560 

585 

608 

1753 

33 

0.95 

401 

0.23 

Net ecosystem exchange 

-156 

-117 

-145 

-419 

29 

0.95 

684 

0.20 

Net primary production 

378 

369 

409 

1156 

20 

N/A 

N/A 

N/A 

Litter production 

235 

242 

265 

742 

51 

N/A 

N/A 

N/A 

Autotrophic respiration 

339 

333 

334 

1016 

18 

N/A 

N/A 

N/A 

Heterotrophic respiration 

221 

252 

264 

737 

36 

N/A 

N/A 

N/A 


Standard deviation on the 3-year analysis in the ensemble of model runs (SD), the fraction of the daily variation of the observations 
explained by the analysis (R 2 ), the number of daily observations (N), and the root-mean-square error (RMSE) of the daily analysis 
over 3 years. 



Total respiration 




Time (days, day 1=1 Jan 2000) 


Fig. 2 Daily observations (o) and analyses (lines) of net 
ecosystem exchange of CO2 (NEE, negative flux = net ecosystem 
uptake), total respiration (R tot ), and gross primary production 
(GPP) for 3 years, 2000-2002. Grey vertical lines indicate the 
standard deviation around the mean of the ensembles, a 
measure of the confidence limits of the analysis. 


and that litter production and heterotrophic respiration 
were of similar magnitude each year. 

The daily analysis of NEE, GPP, and total respiration 
generally closely track their corresponding observa- 
tions (Fig. 2) and have root-mean-square errors (RMSE) 
in the range 0.20-0.32 gCm~ 2 day -1 (Table 4). Most of 
the NEE data lie within one standard deviation of the 
analysis of NEE (Fig. 2, top panel). Similarly, measure- 


3 ~i 



0 365 730 1095 


Time (days, day 1 = 1 Jan 2000) 

Fig. 3 Daily analysis of total heterotrophic respiration (R h , top 
panel), and intermittent observations (symbols, o, with standard 
error bars) and complete analyses (lines) of autotrophic respira- 
tion (Rg, bottom panel) over 3 years, 2000-2002. Grey lines indicate 
the standard deviation around the mean of the ensembles. 


ments of autotrophic respiration all lie within one 
standard deviation of the associated analysis (Fig. 3). 
Cosby et ai (1984) suggest tests for model adequacy 
(i.e., that all the component matrices of the KF, Eqns 
(l)-(4), are correctly specified). The tests analyse the 
innovation sequence specified by 

dj-W) 
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and require that the sequence has a zero mean and 
is serially uncorrelated. Innovations from the model 
forecast of NEE had a zero mean (P > 0.05), but the 
innovations were serially correlated (P<0.01). 

Analysis of C stock dynamics 

The sparse nature of stock data means that the model 
forecast plays a greater role in the analysis, and 
observations correct the model for drift, particularly 
in slow cycling pools. The analysis of foliar C suggests a 
strong seasonal cycle (Fig. 4, top panel), and this is 
partly because of the direct link in the ACM model 
between foliar C (proportional to LAI) and GPP. The 
strong seasonal cycle in observed GPP is transmitted by 
the EnKF to the analysis of foliar C, and the high 
frequency of GPP data generates narrow confidence 
intervals on the analysis of foliar C, compared with 
intervals on the analyses of root or wood C. In three out 
of the four cases, the observed LAI lie within the 
uncertainty of the forecast. The analysis of wood C 
mass indicates a slow gain over 3 years. Each of the 
three observations of wood C lies close to or within one 
standard deviation of the forecast of the previous time 
step (Fig. 4). After an observation is assimilated, the 
uncertainty on the analysis is reduced, and then grows 
slowly with time. The analysis of fine root C suggests 
that mass doubles over the 3 years, and shows how 
reductions in uncertainty when observations are 




Time (days, day 1 = 1 Jan 2000) 

Fig. 4 Daily analysis (lines) and intermittent observations 
(symbols) of selected C stocks. Grey lines indicate one standard 
deviation around the mean of the ensembles. Errors bars (one 
standard error) are also given for each observation. 


assimilated are less than for wood C (Fig. 4), because 
measurement uncertainty on fine roots is larger (Table 
1). For fresh litter, no observations were available, but 
the analysis suggested an annual cycle varying between 
10 and 60gCm~ 2 , with a peak in early spring and a 
minimum in late summer (data not shown). For SOM 
and woody litter, only a single observation was 
available, so the analysis is primarily governed by the 
model forecast (data not shown). The model suggests 
an increase in the SOM/WD over the 3 years. However, 
the limited data mean that confidence limits on the 
analyses for the litter pools are broad. 

Analysis of litterfall 

The analysis suggests that the total production of litter 
over 3 years is 742 gC m~ 2 (Table 4), with fine root litter 
dominating this flux (436 g C m“ 2 ). Because the litterfall 
model is governed by simple turnover rates, and the 
uncertainty on the model was set high, the analyses of 
foliar and woody litterfall are dominated by the 
observations (Fig. 5) and have relatively broad con- 
fidence intervals (Table 4). Lacking any observations. 





Fig. 5 Daily analysis (lines) and intermittent observations 
(symbols) of litterfall. For clarity, errors bars are omitted. 
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the analysis of fine root litter production is driven 
largely by the model, and is, thus, less seasonal. 


Sensitivity analysis 

We explored the sensitivity of the analysis to three 
different factors: (1) the number of time series observa- 
tions assimilated; (2) variation in model error covar- 
iance; and (3) variation in model parameters. 

We explored the sensitivity to the number of 
assimilated time series by reanalysis using the follow- 
ing: (a) the model without any DA, so that model 
forecasts are never corrected with data; (b) assimilating 
only GPP data into the model forecasts; and (c) assi- 
milating GPP data and total respiration data only into 
the forecast (Fig. 6). Without any DA, the model 
forecasts were of NEE were poor (mean of the 
innovations was significantly different from zero, 
P>0.05), and the total NEE prediction over 3 years 
was “251 ±197gCm~ 2 The NEE prediction was a 
clear underestimate of observed summer sink strength. 
The confidence interval of one standard deviation 
encompasses most of the observations (Fig. 6a), 
suggesting that the model error estimate is reasonable. 
Once GPP data were assimilated into the forecast (Fig. 
6b), the analysis was significantly improved. The 
estimate of NEE over 3 years was -413 ± 107 gCnrT 2 , 
and the mean forecast innovations were not signifi- 
cantly different from zero (P<0.05). The analysis 
reproduced the seasonal pattern of NEE, but was 
poorer at predicting the high-frequency variation. 
When both GPP and total respiration data were 
assimilated (Fig. 6c), forecasts of NEE remained 
unbiased (P<0.05), and there was a further reduction 
in the uncertainty on the NEE estimates (over 3 years 
-472 ± 56 g Cm" 2 ). 

We made changes to the values of the model and 
observational error covariances to explore the sensitiv- 
ity of NEE analyses to these parameters (Table 5). A 
halving or a doubling of model error had little impact 
on the analysis. An order of magnitude reduction in 
model error slightly increased sink strength, reduced 
the estimated analysis error, and almost doubled the 
size of the RMSE of analysed NEE vs. observed NEE 
(RMSE). A tripling of model error resulted in a small 
reduction in sink strength and the RMSE. A halving or 
doubling of observational error had only a minor 
impact on NEE predictions and on the comparison of 
analysed with observed NEE. Reducing observational 
error to 10% of the nominal value decreased sink 
strength by only 19gCm~ 2 over 3 years. Tripling 
observational error increased sink strength by just 
6 g C m -2 over 3 years. The sensitivity analyses suggest 



2 H (c) GPP data+fl tot data+ model 



\ 1 i 1 1 1 I ■ 1 l 

0 365 730 1095 

Time (days, day 1 = 1 Jan 2000) 


Fig. 6 Daily analysis (lines) of net ecosystem exchange (NEE) 
generated using (a) model only, no observations; (b) model plus 
gross primary production (GPP) data only; and (c) model plus 
GPP and total respiration data (R to t)- NEE observations are 
shown as open circles. Grey vertical lines indicate the standard 
deviation around the mean of the ensembles. 


that there is little sensitivity in the NEE analysis to 
variations in estimated model and observational errors. 

The nominal parameters and initial conditions used 
in the model (Tables 2 and 3) are not unique, and it is 
important to quantify how uncertainty in these values 
affects the analysis. We used a Monte Carlo technique 
to select alternative parameter sets and initial condi- 
tions, and used these in a multidimensional sensitivity 
analysis. For the model parameters (N = 9) and initial 
conditions (N = 5), we first generated, for each, dis- 
tributions with means equal to the nominal values and 
variances equal to 20% or 50% of the mean. All 
parameters and initial conditions were sampled ran- 
domly from these distributions to generate 400 new, 
unique sets of parameters and initial conditions. We 
generated a new analysis with each of these unique 
parameter sets and initial conditions, and determined 
whether NEE forecasts were unbiased (innovations 
have a zero mean, P>0.05), using the innovation test 
outlined in Results. Of the 400 analyses generated with 
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Table 5 Sensitivity analyses quantifying the impact of 
changes in estimated model and observation error covariances 


Model error 
(% of nominal) 

Observation error 
(% of nominal) 

NEE 

(±SD) 

R 2 

RMSE 

100 

100 

-419 ± 29 

0.95 

0.20 

10 

100 

-440 ±18 

0.81 

0.41 

50 

100 

-423 ± 24 

0.93 

0.25 

150 

100 

-417 ± 33 

0.96 

0.19 

300 

100 

-414 ± 43 

0.96 

0.17 

100 

10 

-400 ± 23 

0.96 

0.17 

100 

50 

-414 ± 26 

0.96 

0.17 

100 

150 

-422 ± 32 

0.94 

0.23 

100 

300 

-425 ± 38 

0.91 

0.29 


Changes were made to either all observational or all model 
errors, using a percentage adjustment to the nominal values 
(see Table 1). The table shows the predicted mean net 
ecosystem exchange (NEE) over 3 years (gCm~ 2 ) ± 1 stan- 
dard deviation of the ensemble. The R 2 and root-mean-square- 
error (RMSE) (gCm -2 ), are for the comparison of the analysis 
with the NEE observations. The first case is the default 
analysis. 

20% variance, 189 had unbiased estimates of NEE, and, 
for the analyses with 50% variance, 75 were valid. For 
the 20% variance analyses, the mean NEE was - 
421 ±17gCm -2 , and for the 50% variance analyses 
the mean NEE was -449 ± 60gCm -2 . Inspection of the 
50% variance cases at the extremes of the range of NEE 
estimation revealed that, while NEE forecasts were 
unbiased, analyses of C stocks (e.g., wood C) were poor. 
The tests of bias and the analyses of NEE were most 
sensitive to changes in the parameter t 2 , the parameter 
that determines autotrophic respiration as a fraction 
of GPP. 

Discussion 

The quality of the analysis 

The tests on the innovations show that the model 
forecast is an unbiased predictor of NEE, although 
autocorrelations in the innovations show that the errors 
are not white. These autocorrelations suggest that Eqns 
(l)-(4) are not perfectly specified, although this might 
be expected given the extreme simplicity of the C 
model we used. Iterative improvements to the model 
should be undertaken in the light of these statistical 
tests of adequacy. For the current analysis, and carefully 
avoiding any prediction into the future, the lack of bias 
in the forecasts means that the results should still 
provide a valid basis for discussion. 

The sensitivity analysis shows how the addition of 
extra data sets, particularly time series of photosynth- 



Fig. 7 Daily variation in the standard deviation around the 
mean of the ensembles on net ecosystem exchange (NEE). 
Symbols show the error on the assimilated estimate, and the line 
indicates the uncertainty assigned to measurement data. 

esis and ecosystem respiration, improves the estimates 
of NEE (compare top panel of Fig. 2 with Fig. 6). The 
confidence interval is approximately halved with 
the assimilation of each additional flux time series. 
With model alone the uncertainty is 196 g Cm -2 ; 
the uncertainty drops to 109 g Cm -2 once GPP data 
are assimilated, to 56gCm -2 with the assimilation 
of respiration data, and to 29 gC m -2 with the assimila- 
tion of NEE data. DA improves on the estimates 
generated by observations alone. The uncertainty on 
daily measurements of NEE was estimated at 
0.5gCm _2 day ' 1 , but the uncertainty on assimilated 
estimates average 0.47 g C m -2 day -1 , and vary depend- 
ing on whether NEE data were available for assimila- 
tion (Fig. 7). Higher uncertainties occurred on days 
when neither NEE nor R tot observations were available 
for assimilation. 

The assimilation process improves the quality of the 
estimates of C dynamics, but the calculated uncertainty 
needs to take account of the full range of uncertainties 
in the model, which can be challenging to quantify. Our 
approach was to use an ensemble of analyses to explore 
how uncertainty in the parameters and initial condi- 
tions of the model affected the quality of the analysis. 
We found that, if all parameters and initial conditions 
were sampled from a normal distribution around their 
nominal values, with a variance set to 20% of the 
nominal values, the mean of the NEE analyses over 
3 years for unbiased models (-421 ± 17gCm -2 ) was 
little different from the nominal analysis 
(-419 ± 29 g Cm -2 ). Increasing the variance to 50% 
broadened the confidence interval to ± 60 g C m -2 , but 
this level of uncertainty in all the model parameters and 
initial conditions is unlikely. Overall, sensitivity analy- 
sis indicated that the uncertainty in model error 
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estimates, and in model parameter estimates, caused 
minor modifications to the analysis. This low sensitivity 
to model errors is important because these errors are 
generally poorly defined. 

For simplicity, we used a daily model and assimilated 
daily data. However, generation of daily data was 
problematical, because flux data are rarely available in 
continuous time series. For NEE data, some hourly gap- 
filling was generally required, and this process neces- 
sarily introduced further uncertainty into the estimate, 
and potentially increased any inherent bias. 

Ecological insights 

The estimates of NEE from our analysis (a net sink) 
differ from those generated for 2001 using more 
conventional approaches for combining flux data 
(Law et al., 2003), which suggested that the forest was 
a small net source of C. NPP estimates from this and 
previous studies are similar, but the analysis of 
heterotrophic respiration for 2001 (184Cm -2 yr -1 ) is 
less than that derived from a more direct estimate based 
on the same automated chamber measurements of soil 
respiration (505Cm” 2 yr _1 ; Irvine & Law, 2002) and 
heterotrophic fractions of soil respiration plus woody 
detritus decomposition (detritus, 67Cm _2 yr _1 ; total, 
301 Cm -2 yr -1 ). This difference is large enough to 
determine whether the forest is a source or sink. 

The strength of the KF approach is that it combines a 
mass balance model with all available data, and thus 
provides a more precise assessment of C exchange than 
pure budget approaches. However, the EnKF relies on 
the model being an unbiased representation of the 
system, and on the data being unbiased. For example, 
the analysis may be degraded by bias in the NEE data, 
because of unreliable night-time measurements. One 
way to remove the night-time bias of eddy flux data 
would be to assimilate hourly data from well-mixed 
periods only. This approach would also avoid the 
problem of using gap-filling to generate daily NEE, but 
would require a model complex model. DA can also 
identify bias, by using the model and other data sources 
to reveal inconsistencies in the analysis. For instance, a 
bias in NEE data would necessarily alter the GPP and 
respiration analyses. Through the interconnections in 
the model, the NEE bias would be transmitted into the 
analyses of GPP and respiration. The transmission of 
bias will thereby increase the differences between the 
analyses of GPP and respiration and the unbiased 
observations of these data. The analysis presented here 
indicates a high degree of consistency among all flux 
data (see R 2 estimates in Table 4), which argues against 
a significant bias. Systematic bias through all observa- 
tions will cause greater problems, as this may conceal 


any inconsistencies. In any DA exercise a careful 
investigation for evidence of bias is vital. 

The model representation of litterfall and turnover is 
oversimplified, and the scarcity of data relating to 
CWD/SOM pools and their turnover means that the 
analysis of these processes is problematical. The 
difference between the EnFK estimate of R h and that 
provided by mass balance suggests that process under- 
standing and measurement capability remain weak in 
this area. Time series data describing the dynamics of 
the soil C and WD pools must be a critical part of any C 
study, so that this area of uncertainty is minimized. The 
broad confidence limits on the modelling of litter 
ensured that the observations determined the values 
of litter fluxes in the analysis. The current model 
formulation is valid for diagnosis, but, until an 
improved litter model is incorporated, it is less suitable 
for prognosis. 

The parameter-fitting exercise suggests an NPP /GPP 
ratio of 53%, slightly larger than the range suggested by 
a survey of multiple forest types (Waring et al, 1998), 
where NPP/ GPP = 0.47 ± 0.04. The parameterized es- 
timates of leaf life span are reasonable (1.0 years), as the 
site is a mix of pine, ferns, and deciduous shrubs. The 
model estimate of fine root life span (1.1 years) is close 
to that suggested by mini-rhizotron data, 1.7 years (Law 
et al, 2001b). However, the turnover rates of woody 
stems (1330 years) and SOM/WD (1033 years) are 
unrealistically large, as a result of the relatively long 
time constants on the dynamics of these pools, 
compared with the length of the observational period 
in this study, and also because of the aggrading nature 
of the forest stand. 

A comparison of NEE estimated by fluxes against the 
overall change in C stocks revealed that the assimilation 
process can break the strict mass balance of the model 
forecast. The analysis suggested that C stocks increased 
by 1990 ± 436 g Cm -2 , a larger, more uncertain sink 
than suggested by summing daily NEE (419 dh 
19 gCm -2 ). The poor match arose because of a gradual, 
but significant accumulation of C in the SOM pool, 
which is poorly constrained and acts as a noise sink. We 
found that increasing the model confidence in the SOM 
prediction improved the mass balance, and reduced the 
accumulation of C in the SOM pool. While the 
mismatch identified here would be problematical over 
very long model runs, in this analysis the SOM pool has 
such a slow turnover that these mass changes have no 
noticeable impact on the predictions of fluxes. For 
prognostic applications, an improved model parame- 
trization is required, and decadal data on woody 
biomass, from inventories and tree ring studies, should 
be assimilated to improve predictions of wood turn- 
over. Field manipulation studies on soil processes will 
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provide insights and constraints into their rates of 
change. 

Operating the model without assimilation, so that 
predictions are uncorrected throughout the 3 years (Fig. 
6a), led to a reduced estimate of the strength of the C 
sink. The model-only simulations revealed a bias in 
NEE forecasts compared with observations. The bias 
largely resulted from the constant, low allocation to 
foliage that leads to a low GPP over the summer 
months. The cumulative effect becomes significant. 
When the model is linked to data via assimilation, the 
GPP and LAI observations ensure a boost of allocation 
to foliage in spring, increased GPP, and a stronger 
summer sink. To improve the prognostic capability of 
the simple model, a detailed leaf phenological model 
with temporal variation in C allocation is required. A 
more effective phenological model should ensure a 
better forecast of maximum LAI in 2001 (Fig. 4). 

The future of DA in ecosystem studies 

Our study has demonstrated why DA provides a 
powerful tool for analysing ecosystem processes, such 
as C cycling. Data alone are often insufficient for this 
task, or at least problematical, because of gaps in time 
series, or methodological uncertainties. Gap-filling 
techniques are generally highly statistical, and fail to 
take account of related data sets or process-based 
models. A conventional modelling approach to ecosys- 
tem analysis can also be problematical. Models are 
often tuned to fit a single data set, and that means the 
model provides little extra information. And if models 
have large numbers of parameters, there may be several 
different combinations of these parameters that can 
be tuned to match the data. Finally, many models 
fail to provide meaningful confidence limits around 
predictions. 

Model DA, as we have demonstrated here with the 
EnKF, solves many of these problems. Assigning 
uncertainties to all observations determines the relative 
importance of data in forming the analysis. All 
predictions have associated confidence limits. The 
EnKF uses a model to generate complete time series, 
so gaps are eliminated. But the gap-filling makes use of 
all other available data at that time, and a modelled 
estimate embodying all prior information from the data 
set that is to be filled. Any increase in uncertainty for 
gap-filled estimates is explicitly calculated. DA makes 
use of all available data sets, so none are wasted in 
model parametrization. In our application here, we 
have estimated the key components of the C cycle from 
a combination of respiration data from chambers, NEE 
data from an eddy flux tower, C mass changes in 
vegetation and soil pools, a detailed model of photo- 


synthesis calibrated from sap flow data, and a mass 
balance model of C dynamics. By attaching error 
estimates to all model components and data sets, the 
KF produces a statistically optimal estimate of NEE 
(and, of course, each state variable in the model). The 
process is an efficient and ordered way of dealing with 
so much data. 

Arguably, the subjective component of the assimila- 
tion process is in choosing the predictive model. By 
keeping this model simple, and by generating the 
model parameters directly from the available data sets, 
with their associated uncertainties, we keep the model 
construction process transparent and reproducible. The 
DA process and the statistical analysis of innovations 
offer the potential to compare alternative models, 
which would improve objectivity. For example, it 
would be intriguing to test a model that predicted 
autotrophic respiration based on the kinetics of meta- 
bolic reaction, rather than indirectly through GPP. 
Examining the analysis to locate when and where 
major corrections occur provides information on model 
reliability. The model we have used here clearly 
oversimplifies seasonal dynamics in litter turnover, 
and foliar and root phenology (Figs 4 and 5), and 
improvements to the model in these areas would be 
worth investigating. Using more detailed models in 
place of the simple model used here should provide 
enhanced analyses. The DA framework provides a 
means of quantifying any such improvement. 

DA offers another important service, the determina- 
tion of minimum data sets required to achieve specified 
confidence levels. The assimilation procedure can 
quantify how changes in the frequency of eddy flux 
data, chamber data, or measures of C pools affect the 
confidence limits on NEE predictions, for example. The 
very dense data collected at the young ponderosa pine 
site are not easily or rapidly reproducible or extensible. 
But it is not currently clear to what degree measure- 
ments could be reduced while still providing a reason- 
able estimate of NEE. DA offers the capability to make 
this determination, and thus play an important role in 
designing efficient regional ecological monitoring sys- 
tems. The results of this initial analysis emphasize the 
importance of time series as constraints. Occasional, 
rare measurements of stocks, such as SOM or CWD 
pools, have limited use in constraining the estimates of 
other components of the C cycle. Long time series are 
particularly crucial for improving the analysis of pools 
with long time constants, such as SOM, woody 
biomass, and WD. Long-running forest stem surveys, 
and tree ring data, offer a rich resource that could be 
assimilated to provide an important constraint on C 
cycling of slow pools. Detailed time series of below- 
ground processes are more challenging to obtain. But 
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both on-line mass spectrometry, to determine continu- 
ously C isotopes of soil C0 2 effluxes, and tunable diode 
laser technology, for monitoring whole ecosystem 
isotopic signatures, may provide powerful sets of 
dynamic data that can be assimilated into models. 
And for extending estimates of NEE across regions, DA 
can play a further important role, by assimilating 
remote-sensing data into the analysis of C cycles. We 
have shown, via sensitivity analysis, how assimilating 
an estimate of photosynthesis - which might be 
provided indirectly by remotely sensed data - improves 
the analysis of NEE. Long-term remote-sensing data 
sets (from the past 30 years) could provide estimates of 
disturbance that could also be assimilated into a model. 
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Appendix 

Carbon dynamic model equations. T is a temperature- 
sensitive rate parameter. 

Fluxes: 


For G, see ACM model below 
= ^2^ 

A( = (1 — h G 

Aw = (1 — f 2 )(l — 13 — U)G 
A r = (1— f 2 )f 4 G 
If = t 5 Cf 
L w — t gC w 
L t = tyCr 

^hi = t s C\i t T 
Rh2 = ^CsomF 
E = fiOitF 

Pools: 


AC* — Af—Lf 
AC W A w L w 
AC r = A t -L t 

AC lit = L f + L w + L r -R hl -D 

ACsom/wd = D - R h 2 

T = 0.5 x exp(0.0693A) (equivalent to a curve with 
Q 10 = 2, T = 1 when A = 10) 
where A is the mean daily air temperature ( °C). 

ACM equations to estimate G 


8c = 


l»d P 

0.5T r + #6^tot ' 


(Al) 


where \fr d is the maximum soil-leaf water potential 
difference (MPa), T r is the daily temperature range ( °C), 
and R to t is the total plant-soil hydraulic resistance 
(MPa m 2 s mmol -1 ). 

V = exp(T max fl 8 ), (A2) 

8c 

where N is the average foliar N (g N m 2 leaf area), L is 
the LAI, estimated by L = Cf/111, T max is the maximum 
daily temperature 


9 =d3 “ fl4, 


(A3) 


2005 Blackwell Publishing Ltd, Global Change Biology, 11, 89-105 


17 


DATA ASSIMILATION OF CARBON DYNAMICS 105 


C- — 1 

'-1—9 
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+ \]{Ca+q ~ p) 2 - 4(C„<j - pa-}) 

a 7 L 2 


eo = 

6 — -0.408 cos | 


L 2 + ag ’ 
f360(D -f 10) re 


365 


180/ 


(A4) 

(A5) 

(A6) 


where <5 is the solar declination (radians) and D the day 
of year. 

s = 24 cos -1 (— tan(A) tan(£))/re, (A 7) 


where s is the day length in hours. If tan(A)tan(£) >1 
then s = 24, A the latitude (radians) 


G = 


g o/j>c(Ca ~~ Q) 
+^c(C« — Cj) 


(fl2S + fl 5 ), 


(A8) 


where G is the GPP in g C m 2 day \ I the irradiance 
(MJm 2 day 1 ). 


Parameters valid for Ponderosa pine GPP predictions 
in ACM model (calibrated from SPA model predictions 
at the young site, Schwarz et al., 2004) 


Parameter 

Value 

a\ 

2.155 

a 2 

0.0142 


217.9 

a 4 

0.980 

as 

0.155 

a 6 

2.653 

a 7 

4.309 

«8 

0.060 

a 9 

1.062 

#10 

0.0006 
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[i] We investigated the relative importance of climatic versus biotic controls on gross 
primary production (GPP) and water vapor fluxes in seasonally drought-affected 
ponderosa pine forests. The study was conducted in young (YS), mature (MS), and 
old stands (OS) over 4 years at the AmeriFlux Metolius sites. Model simulations 
showed that interannual variation of GPP did not follow the same trends as 
precipitation, and effects of climatic variation were smallest at the OS (<10%), largest 
at the MS (>50%), and intermediate at the YS (<20%). In the young, developing 
stand, interannual variation in leaf area has larger effects on fluxes than climate, 
although leaf area is a function of climate in that climate can interact with age-related 
shifts in carbon allocation and affect whole-tree hydraulic conductance. Older forests, 
with well-established root systems, appear to be better buffered from effects of 
seasonal drought and interannual climatic variation. Interannual variation of net 
ecosystem exchange (NEE) was also lowest at the OS, where NEE is controlled more 
by interannual variation of ecosystem respiration, 70% of which is from soil, than by 
the variation of GPP, whereas variation in GPP is the primary reason for interannual 
changes in NEE at the YS and MS. Across spatially heterogeneous landscapes 
with high frequency of younger stands resulting from natural and anthropogenic 
disturbances, interannual climatic variation and change in leaf area are likely to result 
in large interannual variation in GPP and NEE. index terms : 0315 Atmospheric 

Composition and Structure: Biosphere/atmosphere interactions; 1615 Global Change: Biogeochemical 
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1. Introduction 

[2] Temperate forests are an important component of the 
terrestrial carbon cycle, and it is generally believed that 
these forests are net sinks for atmospheric carbon [Myneni et 
al. , 2001 ; Goodale et al, 2002]. Both climate [e.g., Schimel 
et al , 2000, 2001; Schaefer et al, 2002] and time since 
stand-replacing disturbance [e.g., Turner et al , 1995; 
Caspersen et al , 2000; Law et al , 2001b, 2003] have been 
identified as important controls on the net ecosystem 
exchange of C0 2 . Even within a single forest type, dis- 
turbances such as wildfire and timber harvesting have 
resulted in spatially heterogeneous landscapes in the west- 
ern United States [Cohen et al, 1996]. Consequently, many 
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western landscapes are mosaics of forest stands in different 
stages of development (Figure 1) with a range of carbon 
stocks and fluxes [Law et al, 2003, 2004]. 

[ 3 ] Detailed, process-based simulation models provide an 
approach for integrating knowledge of ecophysiological 
processes and scaling these processes to stand and ecosys- 
tem levels [Landsberg and Gower, 1997; Ryan, 2002]. 
Previously, we evaluated the detailed process model, Soil- 
Plant-Atmosphere (SPA), in an old-growth ponderosa pine 
forest (Pinus ponderosa var. ponderosa) in Central Oregon, 
using eddy covariance estimates of whole ecosystem gross 
photosynthesis (GPP) and water vapor exchange [Law et 
al, 2000], and made subsequent modifications of the model 
treatment of root access to soil water [Williams et al, 
2001b]. 

[4] In this study, we took advantage of a unique collection 
of data collected at AmeriFlux sites in young (YS) 20-, 
mature (MS) 90-, and old (OS) >200-year-old stands of 
ponderosa pine in the Metolius River area in central Oregon 
(Figure 1). At these sites, in addition to stand-level data on 
C0 2 and water vapor exchange, the structure, physiological 
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Figure 1. IKONOS image of the Metolius River area of central Oregon showing the locations of the 
young site (YS), mature site (MS), and old site (OS) and the spatial heterogeneity of stand structure and 
development stages. The distribution of stand ages is shown in the lower left. 


characteristics, and soil properties were measured to explain 
variation in ecosystem fluxes in this semi-arid region 
subject to summer drought. 

[ 5 ] The goal of this study was to investigate systemat- 
ically the relative importance of climatic versus biotic 
controls on GPP and water vapor exchange over a 4-year 
period. Our objectives were (1) to parameterize the SPA 
model for three ponderosa pine stands with differing age 
structure and stand characteristics but otherwise growing 
under similar climate and environmental conditions and 
then to compare, for each of the three stands, simulated 
and field-measured estimates of tree transpiration, total 
latent energy (LE) fluxes, and GPP derived from eddy 
covariance data; (2) to cany out a sequence of simula- 
tions designed to evaluate the potential relative impor- 
tance of climate, climatic variation over 4 years, and 
stand-related biotic factors such as leaf area and hydraulic 
conductance as controls on interannual variation of GPP; 
and (3) to examine the implications of interannual vari- 
ation of GPP and ecosystem respiration on net ecosystem 
exchange of C0 2 (NEE). This is the first time we have 
used such a large amount of site-specific data to param- 
eterize the SPA model for different age stands and for 
several years of different climate. We hypothesized that 
interannual variation of GPP at these sites was controlled 
primarily by climatic factors, such as precipitation, but 


that the magnitude of the interannual variation was 
determined by site-specific biotic factors. 

2. Methods 
2.1. Study Sites 

[6] The three sites selected for this study are in the eastern 
Cascade Mountains near Sisters, Oregon, and are located 
within 8 km of one another. The Metolius sites are part of 
the AmeriFlux network [Hargrove et al , 2003]. All three 
sites have an exclusive overstory of ponderosa pine trees but 


Table 1. Selected Site and Stand Characteristics for the Young 
Site (YS), Mature Site (MS), and Old Site (OS), Which are Located 
in the Metolius Area of Central Oregon 



YS 

MS 

OS 

Latitude 

44.437 

44.451 

44.498 

Longitude 

-121.568 

-121.558 

-121,625 

Elevation, m 

1165 

1232 

895 

Mean tree height, m 

4.3 

14.0 

38.9 

Stand age, years 

23 

89 

190 

Foliar nitrogen, g m -2 

2.1 -2.9 

2.1 -2.9 

2.0-3.6 

Sand, % 

54-69 

59-67 

63-67 

Clay, % 

5-11 

8-11 

10-11 

Soil organic matter, % 

2.2 

3.0 

3.0 

Soil porosity, % 

54-59 

53-61 

54-61 
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Figure 2. Monthly climate data for each site during 2002 showing seasonal trends. Vertical bars show 
the range of variation in the monthly values during the 4 years of the study, 1999-2002. 


differ in their age structure and understory composition 
(Table 1). The young pine site (YS) was clear-cut in 1978 
and then allowed to regenerate naturally, and understory 
shrubs account for approximately 40% of the site’s leaf area 
index [Law et al , 2001a]. The mature site (MS) was also 
clear-cut and allowed to regenerate naturally, and the average 
tree age is now about 90 years. The old-growth site (OS) site 
has never been logged and consists of old (250 years), young 
(60-70 years), and mixed-age patches of ponderosa pine. All 
three sites are nearly level, with sandy, well-drained soils 
with relatively little organic matter and little evidence of 
surface runoff (Table 1). Details regarding site conditions and 
soils are described by Law et al. [2001c, 2003]. 

2.2. Climatic Characteristics 

[ 7 ] The three sites have similar climate typical of the semi- 
arid region of central Oregon: warm, dry summers and cool, 
wet winters. Four years of meteorological data (1999-2002) 
were available for the YS and OS, whereas only 1 year of data 
(2002) was available for the MS. The YS and MS are about 
2 km apart and several hundred meters higher in elevation 
than the OS and are ~0.5- 1.0 degrees cooler (Figure 2). At 
the sites, annual precipitation varies between 300 and 
600 mm, with the majority falling between November and 
April as both rain and snow. The winter snowpack generally 


sustains these forests through summer drought, and rain on 
snow events can increase runoff to streams, limiting summer 
water availability to trees. During the 4 years of the study, 
mid-summer rainfall (June 1 to September 30) accounted for 
less than 15% of total annual precipitation. The 30-year mean 
annual precipitation, based on data from the Sisters Ranger 
Station (~15 km east of the study sites), is 360 mm yr“\ 
Although precipitation patterns are highly variable from year 
to year, they are generally synchronous among the three sites. 

Other meteorological conditions are similar among the three 
sites and have substantially less interannual variation than 
precipitation. 

2.3. Environmental Measurements 

[8] In 2002, rates of soil water depletion across different 
soil horizons were estimated from calibrated, periodic time 
domain reflectometry (TDR) measurements in three soil pits 
at each site. Continuous environmental measurements at 
each site included meteorological conditions as well as eddy 
flux measurements of net carbon exchange (F c ) and whole- 
ecosystem water fluxes. Meteorological and micrometeoro- 
logical instruments were installed at a height of 47 m at the 
OS (14 m above the forest canopy), 31 m at the MS (15 m 
above the canopy), and 12 m at the YS (9 m above the 
canopy). Meteorological variables measured at each site “ 
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( consisted of half-hour means calculated from measurements 

at the top of the instrument tower at each of the sites and 
included air temperature, vapor pressure deficit, wind speed 
and direction, global shortwave solar radiation, photosyn- 
, thetically active radiation (PAR), net radiation, and precip- 

itation. After data screening and gap filling, these data were 
I used to drive the SPA model. 

[ 9 ] NEE and total LE fluxes were computed from the 
micrometeorological data using the eddy covariance tech- 
nique. At the YS, eddy covariance measurements have been 
made continuously since the beginning of April 2000 
[Anthoni et al, 2002], and similar measurements have been 
I made at the OS from 1996 until the instruments were moved 

to the MS at the end of 2001 [Anthoni et al , 1999; Law et 
i al, 1999a, 1999b]. Full details of the instrumentation and 

sampling methodology for the OS and YS are described by 
I Anthoni et al [2002] and by M. R. Kurpius et al. (Annual 

carbon exchange along a ponderosa pine chronosequence, 

1 submitted to Journal of Geophysical Research , 2004) (here- 

inafter referred to as Kurpius et al., submitted manuscript, 

) 2004) for the MS. All flux data were summarized and stored 

1 at half-hour intervals, screened, and filtered based on a 

^ friction velocity threshold for wind speed (described below) 

prior to being made available for further analysis [Anthoni 
et al. , 2002]. Data collected on rainy days were excluded 
from the final data set for comparison with the SPA model 
1 because of the effects of rain droplets on the sonic ane- 

mometers and open-path infrared gas analyzers. 

I [ 10 ] Flux-based estimates of gross primary production 

(GPP) were calculated as the sum of NEE and ecosystem 
1 respiration ( R e \ the latter of which was calculated from 

nighttime net carbon exchange from the eddy covariance 
, system and an empirical, Arrhenius-type function of tem- 

perature [Anthoni et al , 2002]. The parameters of the 
1 function were estimated seasonally and thus implicitly 

accounted for gross changes in soil moisture. On the basis 
of examination of F c versus u* at night and comparison of 
R e calculated from nighttime F c (R e e c ) versus R e , calcu- 
lated from summed soil, foliage, and sapwood respiration 
fluxes (R s + R/+ R w ), we determined that local topography 
affects F c under conditions of low turbulence at all three 
sites. We found a dependence ofF c on u* when u* < 0.2 m s~ l 
at YS and OS and u* < 0.35 m s" 1 at MS. Following the 
removal of F c data below the m* threshold and subsequent 
I gap filling using the regression method outlined by Falge et 

al [2001], we achieved good agreement between R c ec and 
R s + Rf+ R w [Anthoni et al , 2002, Kurpius et al., submitted 
manuscript, 2004]. Additional research is being conducted 
r to further understand the role of local topography on 

1 nighttime F c . 

l 

2.4. Structural and Physiological Measurements 

f [ 11 ] Leaf area index (LAI) was estimated at the YS and 

OS in 1999, 2001, and 2002, and at the MS in 2001 and 
* 2002 using the optical approach described by Law et al 

[2001a]. At each site, light transmittance was measured at 
! 39 predetermined locations using an LAI-2000 Plant Can- 

opy Analyzer (Li-Cor, Inc., Lincoln, Nebraska). Mean LAI 
was calculated and then corrected for wood interception and 
for clumping of shoots and needles within shoots according 


to methods outlined by Law et al [2001a]. In 2002, LAI 
was measured in early spring, prior to the expansion of 
current year needles, and again in late summer during 
maximum seasonal leaf area to estimate the seasonal 
changes in leaf area. In 1999 and 2001, LAI was measured 
at the seasonal maximum. Examination of the measure- 
ments of maximum seasonal LAI suggested that leaf area at 
the YS was increasing by about 1 7% per year over the 4-year 
period, whereas maximum seasonal leaf area at the OS and 
MS were essentially constant. 

[ 12 ] Fine root biomass (live roots <2 mm diameter) was 
measured in 2002 to a depth of 1 m by extracting soil cores 
at 18 locations within each site [Law et al, 2003; Sun et al, 

2004]. In the laboratory, roots from the cores were sorted 
into live and dead fractions according to size class. These 
data were used to estimate surface root density (F max ), total 
root biomass (F’totai), and rooting depth (Amx), all required 
parameters for the model (described below). 

[ 13 ] Tree transpiration was estimated by measuring sap 
flux density continuously between April and November at 
each site using a heat dissipation technique [Granier, 1987]. 

Three years of data (2000-2002) were available for the OS 
and YS and one full year of data (2002) was available for 
the MS. Site-level estimates of tree transpiration were 
calculated by scaling the sap flux measurements using 
sapwood area relationships and stand density information. 
Details regarding the methodology and site-level scaling of 
tree transpiration are provided by Irvine et al [2002], 

[ 14 ] Net assimilation ( A-Q curves with an LI-6400; 
LICOR, Lincoln, Nebraska) and foliar nitrogen were mea- 
sured seasonally on shoots from trees and on shrubs. Predawn 
leaf water potential was also measured seasonally, and mid- 
day water potentials were made in mid- summer to estimate 
the critical water potential threshold required by the model. 

2.5. SPA Model 

[ 15 ] The soil-plant-atmosphere (SPA) model is a process 
model that simulates ecosystem photosynthesis, energy 
balance, and water transport. The model was designed for 
direct comparisons with eddy covariance estimates of car- 
bon and water fluxes and has been successfully imple- 
mented in temperate [Williams et al, 1996, 2001b], 
tropical [Williams et al, 1998], and arctic ecosystems 
[Williams et al, 2000]. The model uses 10 canopy layers 
to partition vertical structure of the forest canopy and uses a 
detailed radiative transfer scheme that calculates sunlit and 
shaded fractions of the foliage in each layer [Williams et al, 
2001b]. The model also has 20 soil layers into which fine 
root biomass is distributed that are used to simulate the 
transfer and uptake of soil water received from precipitation 
events. Detailed descriptions of the fundamental equations, 
model logic, algorithms, and the development history of the 
SPA model are given by Williams et al [1996, 2001b], 

2.6. Model Parameterization and Calibration 

[ 16 ] For each site, measured net assimilation in relation to 
leaf internal C0 2 concentrations were used to derive esti- 
mates of maximum carboxylation (F cmax ) and electron trans- 
port (./ m ax) rates [Farquhar and von Caemmerer, 1982], 
which are required parameters of the model (Table 2). The 23 
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Table 2. Key Parameter Values for the Soil-Plant-Atmosphere (SPA) Model Determined From Structural and Physiological 
Measurements at Each of the Study Sites 


Parameter Description 

Units 

YS 

MS 

OS 

Critical leaf water potential prior to cavitation 

MPa 

-1.7 

-1.7 

-1.7 

Daily foliar nitrogen concentration 

g N m 2 

2.1 -2.9 

2.1 -2.9 

1.95-3.6 

Maximum carboxylation capacity (F cmax ) 

pmol m“ 2 s" 1 

48.9 

48.9 

48.9 

Maximum electron transport rate (J m **) 

pmol m~ 2 s -1 

131.6 

131.6 

131.6 

Fcmax per unit leaf nitrogen 

[imol g N -1 s -1 

17.5 

19.6 

14.4 

7 max per unit leaf nitrogen 

pmol g N" 1 s -1 

47.0 

52.6 

38.7 

Daily leaf area index (LAI) in 2002 

-2a 

m m 

0.86-1.46 

2.77-3.44 

2.19-2.74 

Stem-specific hydraulic conductivity 

mmol m 1 s _I MPa -1 

14 

6 

10 

Root resistivity 

MPa s g mmol -1 

12 

3 

10 

Rooting depth 

m 

0 

V© 

1 

© 

1.1 

1.8 

Maximum root biomass per unit volume (F max ) 

g m -3 

3800 

1500 

1400 

Total fine root biomass (F tot aO> min-max 

g m 2 

450-650 

715 

400 


“Denotes m 2 half-surface area of foliage per m 2 ground. 


two optical estimates of LAI during 2002 and measurements 
of foliar nitrogen along with observations of seasonal 
phenology (bud break, needle elongation, and senescence) 
at the sites were used to derive, via linear interpolation, 
daily estimates of LAI and foliar nitrogen for both the 
understory and the trees. Understory leaf area was assigned 
to the lowest canopy layer of the model, and the leaf area of 
the trees was distributed among the upper canopy layers. 
For the YS, additional years of daily LAI and foliar nitrogen 
were prepared in accordance with the observed rate of 
increase of LAI at the site. 

[17] Using field measurements of fine root biomass (F tota i) 
in the top meter of soil, we estimated numerically the 
parameters (F^, D max ) of the exponential function that 
describes the vertical distribution of roots in the model 
[Williams et al , 2001b]. The model also uses depth-specific 
estimates of soil texture to calculate the soil hydraulic 
conductivity and the energy balance at the soil surface. 
These data were available from previous studies for YS and 
OS [Law et al , 2001c] and for the MS [Law et al , 2003]. In 
addition, empirical, site-specific regressions of approximate 
soil water potential as a function of soil water content were 
estimated from TDR measurements and calculated effective 
soil water potential from predawn leaf water potential 
measurements. The remaining model parameters were as 
specified by Williams et al [2001b]. 

[is] Stand hydraulic parameters were calibrated using 
measured estimates of daily tree transpiration derived from 
sap flux density data. We iteratively adjusted the values of the 
two parameters that determine aboveground and below- 
ground hydraulic resistance in the simulated soil-plant- 
atmosphere continuum to determine the combination that 
produced the best agreement between simulated and mea- 
sured tree transpiration. Parameter values were constrained 


based on field measurements of leaf specific conductance 
[Irvine et al , 2004], and final values were selected based on 
the criteria of minimum root mean square error (RMSE) and 
highest coefficient of determination (R 2 ) calculated from 
regressions with measured data. In addition, we calculated 
relative error rates (in comparison with mean measured 
fluxes) and partitioned the mean squared error (MSE) into 
systematic and unsystematic fractions as recommended by 
Wallach and Goffinet [1987, 1989] as cited by Kramer et al 
[ 2002 ]. 

[ 19 ] In summary, sap flux data from each site were used to 
parameterize the SPA model hydraulics, and a single set of 
parameters was determined for each site. Eddy flux data 
were used to evaluate model predictions of total LE and 
GPP. Three years of sap flux data and eddy flux data 
(2000-2002) were available for the YS. For the OS, 3 years 
of sap flux data (2000-2002) were available, but only 
2 years of eddy flux data (2000-2001). For the MS, the 
model was parameterized and evaluated against sap flux and 
eddy flux data for 2002. 

2.7. Simulations 

[ 20 ] To separate relative effects of climatic variation from 
changes in stand structure, for example, LAI and physiology 
(plant water relations), we used the model to perform a series 
of simulations as computational experiments. A base case 
and three experimental scenarios were simulated (Table 3). 
Each scenario was applied to the three sites using the 4 years 
of site-specific climate data (1999-2002), and all 4 years 
were simulated consecutively for a total of 1461 days. These 
4 years represented a 150-200% range of variation of winter 
precipitation (November- April), 300-600% variation of 
summer precipitation (May -October), and over 50% varia- 
tion of annual precipitation (Figure 2). Other climatic vari- 


Table 3. Summary Descriptions of the Four Simulation Scenarios That Investigate Potential Effects of Climatic Variation Versus Stand 
Development and Structure on Forest Carbon and Water Fluxes 3 


Scenario 

Description 

Daily Leaf Area 
Phenology 

Model Parameters 

Soil Moisture 
Limitation 

1 

base case 

site-specific leaf area 

site-specific parameters 

yes 

2 

no drought stress 

site-specific leaf area 

site-specific parameters 

no 

3 

identical LAI and aboveground structure 

MS leaf area 

MS parameters 

yes 

4 

identical aboveground structure and no drought stress 

MS leaf area 

MS parameters 

no 
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, ables showed a range of variation over the 4 years of about 

10-20%. Given the close geographic proximity of the YS 
and MS (~2 km) and the similarity of the climate between the 
two sites, the YS climate data for 1999-2001 (in addition to 
, the MS 2002 data) were used for MS simulations. 

2.7.1. Effects of Climate and Climatic Variation on 
Carbon and Water Fluxes 

[ 21 ] The base case (or first) scenario simulated the effects 
of interannual variation of climate over the 4 years at each of 
the three sites using site-specific LAI, stand structure, and 
hydraulic parameters. The second scenario (no drought 
stress) simulated the removal of soil moisture limitations by 
resetting soil moisture levels to field capacity (~0.3 m 3 m -3 ) 
while using the same leaf area, stand structure, and hydraulic 
parameters as the base case scenario. For the YS, these two 
scenarios were simulated assuming both constant LAI (using 

) the 2002 LAI data) and annually increasing LAI to investi- 

1 gate the added effects of annually increasing leaf area on 

* simulated carbon and water fluxes during the 4 years. Differ- 

ences between the base case and no drought stress scenarios 
i were used to distinguish the effects of soil moisture limita- 

tions (i.e., precipitation) from other climatic effects on carbon 
1 and water fluxes. 

2.7.2. Effects of Aboveground Stand Structure and 

( Canopy Physiology on Fluxes 

[ 22 ] Scenario 3 (identical aboveground stand structure 
and physiology) simulated the effects of replacing the entire 
aboveground structure of the YS and OS with the above- 

I ground structure (LAI, tree heights, vertical distribution of 

LAI) and physiology (canopy hydraulic conductance) of the 
' MS. Differences between the base case (scenario 1) and 

scenario 3 were used to help determine the effects of 
aboveground stand structure and canopy physiology on 
carbon and water fluxes. Variations of this scenario were 
used to distinguish the effects of LAI from other compo- 
nents of stand structure and physiology. 

2.7.3. Effects of Root Density and Rooting Depth on 
Fluxes 

[ 23 ] Scenario 4 (identical aboveground stand structure 
and no drought stress) was identical to scenario 3 except 

| that soil moisture limitations were removed (like scenario 2). 

Under this scenario, virtually unlimited access to soil water 
effectively minimizes differences in the rooting character- 
istics among the sites. Thus, differences between scenario 3 
1 and scenario 4 were used to assess effects of differences in 

1 soil properties and differences in root density and distribu- 

1 tion among the sites. 

[ 24 ] The scenarios are summarized in Table 3. For sce- 

t narios 2-4, the MS was selected as the reference site 

because it represents the upper limit of LAI and productivity 
among stands in the region [Law et al., 2003]. Consequently, 
comparisons with the MS under these scenarios show the 
I maximum potential effect of these simulated changes in 

stand structure on carbon uptake and water vapor exchange. 

2.8. Annual Estimates of NEE 

1 [ 25 ] Annual estimates of NEE were calculated as the 

difference between total annual GPP (from the simulations) 
and annual estimates of R e . Estimates of R e for each site 
were calculated as the sum of R s , determined from auto- 


mated chamber measurements [Irvine and Law , 2002], and 
empirically modeled estimates of Rf and R w [Law et al ., 

1 999b], plus site-specific estimates of mean annual fine and 
coarse woody debris decomposition [Law et al., 2003]. Four 
years of sapwood volume data were used to calculate R w at 
the YS. In the absence of similar data for the MS and OS, 

R w was assumed constant over the 4 years of the study. 
Previous studies at the YS and OS indicate that R w is small 
fraction of R e (<6%) [Law et ah, 1999b, 2001c], and this 
assumption had little impact on the annual estimates of R e . 

3. Results and Discussion 

3.1. Model Parameterization 

[ 26 ] We were able to identify a single pair of hydraulic 
parameters for each site that was held constant during the 
multiyear simulations. In contrast, site -specific LAI and 
total foliar nitrogen varied seasonally. On the basis of the 
2002 measurements, total (overstory plus understory) LAI 
ranged from 0.86 to 1.46 nr half-surface area of foliage per 
m 2 ground at the YS, 2.19 to 2.74 m 2 m~ 2 at the OS, and 
2.77 to 3.44 m 2 m~ 2 at the MS (Table 2). The parameter- 
izations indicated that rooting was deepest at the OS (1 .8 m) 
and shallowest at the YS (1.0 m). However, fine root 
biomass in the top meter of soil was highest at the MS 
and lowest at the OS. Rooting depth and root biomass were 
assumed constant at the MS and OS during the simulations, 
but at the YS, because of its stage of development, fine root 
biomass increased over the 4 years from 450 to 650 g m -2 , 
and rooting depth increased from 0.9 to 1.0 m. In terms of 
water transport, the most conductive stems were at the YS, 
and the least conductive stems were at the MS. The most 
resistive roots were at the YS, and the least resistive were at 
the MS. Leaf specific conductance, which combines the 
stem and root hydraulic parameters, at the top of the forest 
canopy was 1.70, 0.28, and 0.20 mmol s” 1 m -2 MPa 1 at 
the YS, MS, and OS, respectively. 

3.2. Model Comparison With Sap Flux Measurements 

[ 27 ] Simulations of the YS, MS, and OS indicate that the 
SPA model can accurately predict seasonal patterns of tree 
transpiration estimated from measured sap flux using simple 
parameterizations (Table 4). The SPA model explained 76% 
to 87% of the daily variation of tree transpiration for years 
2000 through 2002 across the three sites. RMSE varied little 
(0.1 1 -0.20 mm d" 1 ) among the seven site-year simulations, 
suggesting that it may be difficult to resolve model-data 
discrepancies of less than about 0.1 mm d -1 . The average 
relative error (calculated as RMSE/mean daily measured sap 
flux) between simulated and measured transpiration was 24% 
(±0.14 mm d -1 ). Furthermore, after partitioning the mean 
square error into systematic and unsystematic fractions, on 
average more than half (57%) of the model-measurement 
discrepancies can be attributed to unsystematic (or random) 
errors, suggesting that the model is accurately capturing the 
daily dynamics of tree transpiration (Figure 3). 

3.3. Model Comparison With Latent Energy Fluxes 

[ 28 ] Comparisons between simulations of total daily 
latent energy (LE) fluxes at each of the three sites and 25 
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Table 4. Comparisons of Daily Estimates of Simulated Tree Transpiration, Total Latent Energy Exchange (Total LE), and Gross Primary 
Production (GPP) With Measured Values Derived From Sap Flow and Eddy Covariance Measurements 8 

Tree Transpiration Total LE GPP 

Rel. Rel. Rel. 


Site 

Year 

R 2 

RMSE, 
mm d 1 

Error, 

% 

Bias, 

% 

R 2 

RMSE, 
mm d -1 

Error, 

% 

Bias, 

% 

R 2 

RMSE, 
g C m -2 d -1 

Error, 

% 

Bias, 

% 

YS 

2000 

0.82 

0.14 

28 

17 

0.79 

0.36 

30 

42 

0.90 

0.73 

24 

85 


2001 

0.87 

0.13 

27 

13 

0.81 

0.32 

32 

17 

0.86 

0.56 

23 

48 


2002 

0.86 

0.20 

30 

48 

0.80 

0.38 

33 

24 

0.90 

0.49 

17 

19 


Avg 

0.85 

0.16 

28 

26 

0.80 

0.35 

32 

28 

0.89 

0.59 

21 

51 

MS 

2002 

0.84 

0.15 

19 

32 

0.50 

0.62 

50 

67 

0.46 

0.74 

12 

20 

OS 

2000 

0.85 

0.11 

19 

54 

0.62 

0.51 

36 

81 

0.40 

0.56 

13 

38 


2001 

0.87 

0.11 

23 

65 

0.64 

0.36 

31 

49 

0.67 

1.13 

22 

91 


2002 

0.76 

0.14 

22 

67 

- 

- 

- 

. 

- 

- 

_ 



Avg 

0.83 

0.12 

21 

62 

0.63 

0.44 

34 

65 

0.53 

0.84 

18 

64 

Overall Avg 


0.84 

0.14 

24 

43 

0.70 

0.42 

35 

47 

0.70 

0.70 

19 

52 


a RMSE is root mean square error, and the relative error (Rel. Error) is the RMSE as a percentage of the daily mean value determined from field 
measurements. Bias is the fraction of the mean square error that can be attributed to systematic (as opposed to random) deviations from the measured 
values. Dashes indicate no data available. 


eddy covariance measurements indicate that the model 
was able to explain 50% (MS) to 81% (YS) of the daily 
I variation of total measured LE among the six site-year 

\ combinations with eddy flux data (Table 4). The average 

relative error between the model and the flux measure- 
ments was about 35%, and the overall average RMSE 
was ±0.42 mm d -1 . Relative error rates as well as the 
I proportion of the error attributed to systematic differences 

were smallest at the YS and higher at the MS and OS. 
Eddy covariance-based estimates of integrated total LE at 
| the YS and MS during non-rain days were higher than 

\ simulated estimates (Table 5), although the simulated 

estimates were within the uncertainty limits for the YS. 
Because tree transpiration was the largest component of 
total LE (>50% at each site), the higher relative error 
rates compared with tree transpiration were most likely 


associated with other components of total LE such 
as shrub transpiration, wet canopy evaporation, or soil 
evaporation. 

[ 29 ] Because the same model hydraulic parameters were 
applied to both the understory and overstory canopy layers, 
we would have expected larger discrepancies between the 
model and field measurements at the YS (understory LAI 
40% of total). However, relative error rates (or discrepan- 
cies) for total LE at the YS were smaller than the 
corresponding relative error rates at either the OS or MS, 
contrary to our expectations (Table 4). Thus it is unlikely 
that shrub transpiration, alone, could explain the discrep- 
ancies in simulated versus measured total LE, especially at 
the OS and MS. 

[ 30 ] Evaporation from soil surfaces could also be con- 
tributing to the discrepancies between the flux measure- 


Modeled versus measured tree transpiration (2002) 



Figure 3. Simulated versus measured tree transpiration during 2002, the only year with measured sap 
flux data for each site. The 1:1 line indicates that the simulated values are relatively unbiased. 


26 
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Table 5. Integrated Water Vapor Fluxes and Annual Carbon Fluxes Derived Eddy Covariance Measurements and Simulations* 


Site (Year) 

LE, b 

mm 

GPP (g C m~ 2 yr 1 ) 

NEE, C g C m~ 2 yr' 1 

Measured 

Modeled 

Measured 

Modeled 

Measured 

Modeled 

YS (2001) 

224 ± 36 

199 

704 ± 168 

674 

45 ±6 

48 

YS (2002) 

269 ± 43 

234 

809 ± 193 

828 

113 ± 15 

192 

MS (2002) 

296 ±51 

241 

1360 ft 333 

1539 

371 ±53 

516 


“For the OS, insufficient eddy flux data were available for 2001 to estimate annual fluxes, and no eddy flux data were available for 2002. 
b Rain days were excluded from the integrated estimates of LE because of the effects of rain droplets on the eddy covariance instrumentation. Thus the 
LE estimates are not true annual estimates. 

c Modeled estimates of NEE are determined from simulated estimates of GPP plus empirical estimates of ecosystem respiration derived from 
measurements of soil CO 2 efflux, foliage respiration, sapwood respiration, and woody litter decomposition (see section 2.8). 


ments and the model. Williams et al [2001b] reported 
good agreement between the SPA model and measure- 
ments of soil evaporation at the OS during August 1997. 
In the current study, however, we utilized field measure- 
ments of tree transpiration and expanded the number of 
sites and the number of years of data used for the model- 
measurement comparison. The largest discrepancies be- 
tween the eddy covariance-based estimates of total LE 
and the model occurred usually during June and July, 
when the model underestimated total LE by ~1 mm d _f 
or >40% at the MS. Simulated soil evaporation at the MS 
during the summer months was relatively constant and 
averaged ~0.2-0.3 mm d _1 . Although simulated soil 
evaporation compared favorably with LE fluxes measured 
by a subcanopy eddy covariance system installed at the 
MS during August- September 2002, no measurements 
were available from earlier in the summer when measure- 
ment-model discrepancies were larger. Thus, to further 
improve the model’s capability to simulate soil evaporation 
at these sites, additional measurements of soil evaporation 
will be needed that span the growing season. 

[ 31 ] Another possible source for discrepancies between 
simulated and measured total LE is differences in the 
sampling footprint used to calculate tree transpiration from 
sap flow measurements and the sampling footprint for the 
eddy covariance measurements, even though a rigorous 
approach was used to scale the sap flow measurements to 
site [Irvine et al., 2002]. Scaled estimates of transpiration 
are sensitive to the measured basal area and sapwood area 
relationships in the field plots, and although the field plots 
were chosen to be representative of each site, structural 
characteristics of the plots may not match those of the area 
being integrated by the eddy covariance technique. Other 
recent studies have also reported difficulties in reconciling 
sap flow measurements with latent energy flux measure- 
ments using the eddy covariance technique [Wilson et al, 
2001; Meiresonne et al, 2003], and we submit that the 
footprint issue requires further investigation. 

3.4. Model Comparison With GPP Estimates Derived 
From Eddy Flux Data 

[ 32 ] The average relative error between modeled and 
measured estimates of daily GPP among the three sites 
was 19%, and RMSE varied by about a factor of 2 from 
0.49 g C m -2 cT l to 1.13 g C m -2 d _1 (Table 4). 
Simulations of the YS, though, explained considerably more 
day-to-day variation in GPP than at either the MS or OS, 
and roughly half of the mean square error could be 


attributed to systematic differences between the modeled 
and measured values. Although the overall relative error 
between the eddy covariance-based estimates and simulated 
estimates of GPP was generally good (19%, Table 4), the 
lower R 2 values and higher bias indicate that the model 
misses some of the day-to-day variation associated with the 
carbon dynamics of the sites. In general, though, as time 
series data are aggregated over longer periods, discrepancies 
between predicted and measured values are substantially 
reduced because differences due to unsystematic (or ran- 
dom) variation tend to cancel one another [Goulden et al, 
1996; Moncrieff et al, 1996]. 

[ 33 ] Annual estimates of simulated GPP were well within 
the uncertainty limits of the flux measurements (Table 5). 
Simulated estimates were higher than the measured esti- 
mates in 2002 but lower than the measured estimate for the 
YS in 2001. Thus, although the model may be under- 
estimating total LE at the study sites, the model appears 
suitable for simulating seasonal and annual responses of 
transpiration and GPP to changes in climate and climatic 
variation. 

3.5. Simulation Scenarios 

3.5.1. Effects of Climate and Climatic Variation 

[ 34 ] Precipitation at the sites varied ~55% over the 
4 years: 1999 was the wettest year followed by 2001, 

2000, and 2002 was the driest year. When LAI at the YS 
was held constant over the 4 years to separate effects of 
climatic variation from changes in leaf area, GPP was 
highest during 1999 and lowest in 2002, following interan- 
nual trends in total annual precipitation. At the MS and OS, 
however, annual GPP was highest in 1999 and lowest in 

2001, which mirrored the interannual pattern of soil water 
recharge from winter precipitation. 

[ 35 ] Despite similar climate across the three sites, the base 
case simulations (scenario 1) showed that interannual var- 
iation of GPP differed by site in response to climatic 
variation over the 4-year period (Table 6, Figure 4, base 
case). Annual GPP in the base case simulations varied 
<10% per year at the OS to >50% per year at the MS 
across the 4 years, implying that interannual variation of 
GPP among the 4 years was highest at the MS and lowest at 
the OS. Under the base case scenario with constant LAI, 
interannual variation of GPP at the YS across the 4 years 
was 19%. The interannual variation of GPP at the MS was 
similar to the interannual variation of precipitation, and the 
interannual variation of GPP at the YS was comparable to 
that of the MS when simulations increased the LAI at the 
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Table 6. Annual Summaries of the Simulated Carbon and Water Vapor Fluxes From the Simulations 8 


Scenario 


Total LE, mm yr 1 

Tree Transpiration, 
mm yr -1 

GPP, gCm‘ 2 y 

-1 

YS 

MS 

OS 

YS 

MS 

OS 

YS 

MS 

OS 

Base case 

mean 

309 | 279 

304 

299 

136 | 115 

171 

161 

801 | 679 

1401 

1098 


std dev 

31 | 17 

37 

13 

13 | 16 

34 

9 

55 | 93 

219 

30 


percent variation 

33 | 17 

40 

13 

28 | 38 

77 

17 

19 | 45 

53 

7 

No drought stress 

mean 

434 | 385 

442 

327 

241 | 204 

309 

187 

1013 | 846 

1886 

1168 


std dev 

18 | 32 

20 

14 

4 | 32 

11 

6 

8 | 132 

49 

15 


percent variation 

11 | 17 

12 

11 

5 | 52 

10 

9 

2 | 55 

7 

3 


percent change 

40 | -11 

45 

9 

78 | -15 

80 

16 

26 | -17 

35 

6 

Identical LAI and aboveground structure 

mean 

305 

304 

369 

191 

171 

227 

1413 

1401 

1648 


std dev 

25 

37 

20 

23 

34 

22 

100 

219 

83 


percent variation 

23 

40 

16 

37 

77 

31 

21 

53 

14 


percent change 

-1 

- 

23 

40 

- 

42 

76 

- 

50 

Identical aboveground structure and no drought stress 

mean 

392 

442 

438 

283 

309 

289 

1686 

1886 

1816 


std dev 

20 

20 

10 

9 

11 

10 

24 

49 

30 


percent variation 

14 

12 

5 

9 

10 

10 

4 

7 

4 


percent change 

27 

45 

46 

108 

80 

27 

110 

35 

65 


Percent variation describes the relative range of variation over the four simulated years, 1999 -2002. 

“Percent change describes the relative difference between a given scenario and the base case scenario. Dashes indicate not applicable. The vertical bars 
for the YS separate the results for the constant LAI simulations on the left from the annually adjusted LAI simulations on the right. 


YS to match the measured LAI. In contrast, the interannual 
variation of GPP at the OS was much less than the 
interannual variation of precipitation. Interchanging climate 
driver data among the sites had little effect on mean annual 
GPP at either the OS or MS and reduced mean GPP at the 
YS by no more than 5%. 

[36] Seasonally, interannual variation of GPP at the MS 
was most pronounced during July -September, while varia- 
tion at the YS was most noticeable earlier during May -July 
(Figure 4, base case). At the OS, the interannual variation 
appeared uniformly distributed throughout the year. 

[ 37 ] Base case simulations (scenario 1) of tree transpira- 
tion and total LE showed noticeably more interannual 
variation than GPP because of larger interannual variation 
in LE and precipitation earlier in the year (Table 6, Figure 5, 
base case). Among-site differences in transpiration were 
substantially smaller than the corresponding differences in 
GPP (Table 6): about 25% for mean annual tree transpira- 
tion (assuming constant LAI at the YS) and less than 5% for 
mean annual total LE. 

[ 38 ] Increased soil water availability at the MS and YS, 
under scenario 2 (no drought stress), resulted in large 
increases in mean annual tree transpiration of >75% and 
increased total annual LE by >40% (Table 6). Moreover, 
the increased tree transpiration and total LE were most 
evident during the period of otherwise maximum drought 
stress (Figure 5). The simulations suggest that increased soil 
moisture from precipitation during the growing season 
would probably have a greater effect on carbon and water 
fluxes in these forests than increased precipitation at other 
times of the year [Goldstein et al , 2000; Xu and Baldocchi, 
2004]. Mean simulated tree transpiration and total LE at the 
OS, on the other hand, increased only by about 16% and 
9%, respectively, reflecting a smaller effect of increased soil 
moisture in the upper soil layers at this site. These results 
are consistent with those of Irvine et al. [2004], who 
determined that approximately 79% of water used during 
the summer months of 2002 at the YS was extracted from a 
depth of 80 cm or less, whereas almost half (47%) of the 


water extracted at the OS during the same months came 
from below 80 cm depth. 

[ 39 ] Simulations that increased soil water availability and 
eliminated seasonal drought stress (scenario 2) increased 
mean annual GPP at the YS and MS by 26% and 35%, 
respectively (Table 6). Furthermore, simulations clearly 
show that the increased GPP occurs during the late summer 
and early fall months when drought stress would normally be 
experienced (Figure 4). In contrast, the simulated response of 
the OS was small: Mean annual GPP increased by 6%, and 
the seasonal pattern of GPP was similar to that of the base 
case scenario. Therefore the benefits of increased soil water 
availability had a much bigger effect at the Y S and MS than at 
the OS. 

[40] Increased soil water availability subsequently re- 
duced interannual variation of GPP to similarly low levels 
at all three sites (2-7%, Table 6). Hence the simulations 
confirm that interannual variation of precipitation is an 
important causal factor controlling GPP and probably 
accounts for most of the interannual variability of GPP in 
these forests. Under scenario 2, though, after removing the 
variation of rainfall, the dominant sources of interannual 
climatic variation were radiation inputs and temperature, 
and the interannual variation of radiation probably had a 
larger effect on the interannual variation of GPP at the MS 
(7%) because it had highest LAI among the three sites. In 
contrast, the YS experienced the least interannual variation 
of GPP (2%, assuming constant LAI), and although the 
interannual variation of solar radiation was higher at the YS 
compared to the OS, even during the year with the lowest 
solar radiation inputs (2001), there was sufficient radiation 
to effectively saturate the lower LAI canopy of the YS. 

[41] Other studies reached similar conclusions regarding 
effects of interannual variation of climate on carbon uptake 
in forests. For example, Barford et al [2001], in a 9-year 
study of eddy covariance data from a deciduous hardwoods 
stand in Massachusetts, attributed interannual variation of 
carbon uptake and net carbon exchange to the variation of 
climatic characteristics such as cloudiness and low temper- “ 
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month 


MS 



1 2 3 4 5 6 7 8 9 10 11 12 

month 


OS 



1 2 3 4 5 6 7 8 9 10 11 12 

month 


Figure 4. Monthly GPP under the base case and no drought stress scenarios showing effects of seasonal 
drought at the sites. Differences between the two curves are an estimate of the reduction in GPP caused 
by drought. Vertical bars show the range of variation over the four simulated years, 1999-2002. 

atures, which influenced growing-season length. Griffis et [ 43 ] Increasing the LAI at the YS to match interannual 
al. [2003], based on 6 years of eddy covariance data, changes in measured values during the base case simula- 
associated higher GPP with warmer climatic conditions at tions produced higher interannual variation of the fluxes of 
an aspen site in Saskatchewan, Canada. However, to our carbon and water vapor (Table 6). For GPP, interannual 
knowledge, our study is the first to suggest that gross carbon variation increased from 19% to 45%, and the pattern of 
uptake of nearby stands of the same forest type, with largely GPP relative to the pattern of annual precipitation was 
similar climate, can have substantially different sensitivities reversed. GPP was highest in 2002 and lowest in 1999, 
to interannual variation of climate. which reflected the importance of stand leaf area in deter- 

3.5.2. Effects of Aboveground Stand Structure and mining gross carbon uptake. Similar responses to increasing 
Canopy Physiology LAI during the 4-year simulations were also noted for tree 

[ 42 ] Mean annual simulated GPP during 1999-2002 was transpiration and total LE. Under scenario 2 (no drought 
highest at the MS and lowest at the YS, consistent with the stress), the interannual variation of GPP over the 4-year 
rank differences in LAI among the sites (Table 2). Mean period increased despite the stabilizing effect of increased 
annual GPP was 1401, 1098, and 679 g C m -2 yr” 1 at the soil water availability (Table 6). Thus the simulations 
MS, OS, and YS, respectively (Table 6, base case). When suggest that increases in leaf area, alone, may have a 
LAI was held constant at the YS (using the 2002 LAI), tendency to increase the sensitivity of a stand to interannual 
mean annual GPP increased by ~18% to 801 g C m' 2 yr" 1 . variation of climate, at least in low LAI systems like our ^ 
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YS 



MS 



OS 



Figure 5. Monthly tree transpiration under the base case and no drought stress scenarios, showing 
effects of seasonal drought at the sites. Differences between the two curves are an estimate of the 
reduction in GPP caused by drought. Vertical bars show the range of variation over the four simulated 
years, 1999-2002. 


sites. In stands with higher LAI, though, changes in leaf 
area could produce feedbacks in other components of the 
stand’s overall energy balance. 

[44] Scenario 3 (identical aboveground structure and 
physiology) replaced the aboveground structure of the YS 
and OS with the LAI, foliage distribution, and canopy 
hydraulic conductance of the MS. Most of among-site 
differences in GPP can be explained by differences in stand 
structure and canopy physiology (Figure 6, third plot). 
Compared to the base case, mean annual GPP increased 
76% and 50%, respectively, at the YS and OS, and differ- 
ences in mean annual GPP among the three sites were 
reduced to 1 8% (Table 6). Additional simulations at the OS, 
which separated the effects of LAI, hydraulic conductance, 
and vertical canopy structure, revealed that 16% of the 50% 
increase in mean GPP could be attributed to the increased 


LAI, and most of the remainder was caused more by the 
simulated changes in height and vertical distribution of 
foliage rather than the simulated changes in hydraulic 
properties of the sites. In contrast, at the YS, 61% of the 
76% increase in mean GPP was related to the increased 
LAI, with most of the remainder being attributed the 
simulated changes in height and vertical distribution of 
foliage, which increased absorption of radiation by the plant 
canopies. 

[ 45 ] Simulated responses of tree transpiration and total 
LE to changes in aboveground stand structure under 
scenario 3 differed somewhat from those of GPP (Table 6, 
Figure 7). Mean annual tree transpiration at the YS and 
OS increased by ~40% relative to the base case scenario, 
but in contrast with simulated GPP, differences among the 
sites in mean annual transpiration increased as well, suggest- ^ 
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1 . base case 2. no drought stress 



3. identical aboveground structure & physiology 



4. identical aboveground structure, physiology 8> no drought stress 



Figure 6. Comparisons of monthly GPP among the four simulated scenarios for each of the three sites. 
Vertical bars show the range of variation over the four simulated years, 1999-2002. The simulations for 
the YS assume that LAI was constant over the 4 years. 


ing that belowground differences were contributing among- 
site differences. 

[ 46 ] Most of the seasonal variation of the fluxes of carbon 
and water vapor was evident during the months of July- 
October (Figures 6 and 7), similar to the base case scenario. 
Furthermore, despite identical aboveground structure at 
each of the sites, there were still distinct differences in the 
interannual variation of GPP among the sites, which were 
53%, 21%, and 14% at the MS, YS, and OS, respectively 
(Table 6), further suggesting that belowground site differ- 
ences were responsible for the pattern of interannual vari- 
ation. Similarly, the same pattern of interannual variation 
among the three sites was evident for tree transpiration and 
total LE. 

[ 47 ] Other simulated changes in aboveground structure 
besides LAI influenced the hydraulic properties of both 
stands but had a bigger relative impact on GPP at the OS 
than at the YS. The net effect of altering the tree heights 
and hydraulic conductance increased overall canopy hy- 
draulic conductance at the OS but decreased canopy 
conductance at the YS. Hydraulic limitations have been 
proposed as an explanation for reduced productivity in 
older trees [e.g., Yoder et al , 1994; Ryan and Yoder , 
1997], but see Ryan et al [2004]. The effects of these 
simulated changes are a logical extension of previous re- 
search on hydraulic limitations [e.g.. Waring and Running , 
1978; Zimmermann, 1983; Yoder et al, 1994; Ryan and 


Yoder ; 1997] and have been discussed previously by Williams 
et al [2001a]. 

[48] Our simulations also suggest that rapidly developing 
stands, such as the YS, may be prone to higher interannual 
variation of GPP than stands with stable LAI (Table 6, base 
case). Moreover, the interannual variability of GPP in- 
creased in the simulations that removed soil moisture 
limitations at the YS (no drought stress scenario). In the 
absence of water limitations, GPP is constrained on an 
annual basis mostly by the amount of PAR intercepted by 
the foliage, and over the 4-year simulation, the fraction of 
absorbed PAR at the YS increased by about 40%, which 
explains most of the increased variability of GPP despite the 
otherwise stabilizing effect of increased soil water avail- 
ability. After accounting for the increase in LAI, the 
remaining variation (15%) was similar to the variation of 
the base case scenario with stable LAI (19%) and can be 
attributed to interannual climatic variation. In general, 
though, our results suggest that following disturbance, 
young, rapidly developing stands may be more sensitive 
to variations in climatic conditions than stands with more 
stable LAI, which, in turn, could cause these stands to shift, 
on an annual basis, from net sources or net sinks of carbon, 
assuming that the interannual variation of GPP is larger than 
the interannual variation of ecosystem respiration [Law et 
al, 2001c, 2003], This possibility is discussed in more 
detail below. 
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Figure 7. Comparisons of monthly tree transpiration among the four simulated scenarios for each of the 
three sites. Vertical bars show the range of variation over the four simulated years, 1999-2002. The 
simulations for the YS assume that LAI was constant over the 4 years. 


[ 49 ] Our purpose for simulating changes in stand struc- 
ture, though, was to investigate possible explanations for the 
differential response and sensitivity among the sites to 
increased precipitation and interannual climatic variation, 
and although the simulated changes in aboveground stand 
structure and physiology had a large impact on GPP at the 
sites, they had little effect on the interannual variation of 
GPP. The range of interannual variation among the sites 
(with identical aboveground structure) was similar to the 
range of variation under the base case scenario, and the rank 
differences among the three sites were the same (Figure 6, 
Table 6). Hence we concluded that differences in above- 
ground stand structure and physiology, alone, could not 
explain the site- specific differences in sensitivity to inter- 
annual climatic variation. 

3.5.3. Effects of Root Density and Rooting Depth 

[ 50 ] Under scenario 3 (identical aboveground structure 
and physiology), the only differences between the sites were 
those differences related to belowground structure (fine root 
density and distribution) and to physical soil properties, 
which were minor (Table 1), and the simulations under 
scenario 4 (identical aboveground structure and no drought 
stress) effectively minimized differences in carbon uptake 
and water vapor exchange among the three sites (Figure 6). 
Virtually unlimited soil water availability increased mean 
annual GPP at the YS and OS by 19% and 10%, respec- 
tively, compared with the previous scenario, and the vari- 


ation of mean annual GPP among the three sites was 
reduced to 12% (Table 6). Furthermore, the interannual 
variation of GPP was similar among the sites and ranged 
from 4 to 7%. Hence the simulated increase in soil water 
availability compensated for differences in belowground 
structure among the sites and effectively eliminated differ- 
ences in root density and distribution. Additional simula- 
tions using identical climate data for each site (rather than 
site-specific data) produced virtually identical interannual 
variation, 7.5 -7.7% among the three sites. Thus most of the 
differences in interannual variation of GPP among the sites, 
evident in the previous simulations, could be attributed to 
site-specific differences in fine root density and root distri- 
bution. Therefore the smaller interannual variation of GPP 
at the OS evident in the base case simulations (Figure 6) can 
be explained by root systems at the site that can access 
water deeper in the soil profile. 

[ 51 ] The effects of scenario 4 on tree transpiration and 
total LE were similar to those of GPP. Mean annual fluxes 
increased and differences among the sites decreased 
(Figure 7). Differences in interannual variation of transpi- 
ration and total LE among the sites were also reduced 
(Table 6), and although most of the among-site differences 
of total LE were related to differences in tree transpiration, 
some of the differences were also related to soil evaporation. 

[ 52 ] An alternative explanation for lower interannual 
variation of GPP at the OS and the site’s modest response 3. 2 
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to simulated increased soil water availability could be that 
soil moisture levels at the OS were higher than at either the 
YS or MS. Precipitation, though, at the three sites was 
similar (Figure 2), and measurements of soil water content 
during the growing season in the topmost meter of soil 
indicated that soil moisture was also similar among the three 
sites (J. Irvine, unpublished data, 2002). Furthermore, 
measured and simulated soil water content values were also 
comparable. However, predawn foliage water potential 
measurements during the 1999-2002 growing seasons 
showed that the OS had higher levels of available soil water 
than either the YS or MS, which suggested that OS 
experienced less drought stress than the other two sites 
[Irvine et ai, 2004]. Consequently, among-site differences 
in soil characteristics and soil water content in the top meter 
of soil were minor and do not explain the sites’ differential 
response to interannual climatic variation or a simulated 
increase in precipitation. 

[ 53 ] Simulation scenarios 3 and 4 provide compelling 
evidence that the most likely explanation for both the 
modest reduction in GPP at the OS in response to drought 
as well as the site’s lower interannual variation of GPP is 
that tree roots at the site are able to access soil water from 
deeper in the profile, whereas tree roots at the YS and MS 
cannot. In addition to the simulation results, there are some 
empirical data that support the conclusion of a greater 
rooting depth at the OS. Irvine et ai [2002] concluded that 
a deeper rooting depth at the OS (compared to the YS) was 
the most likely explanation for less negative soil water 
potentials at the OS in an empirical study of sap flux 
measurements and soil C0 2 fluxes at the YS and OS. 
Furthermore, analyses of stable oxygen isotope data, al- 
though not conclusive, indicated that the 6 18 0 signature of 
xylem water in trees at the OS was more similar to that of 
nearby spring water than the corresponding isotope signa- 
ture of soil water in the surface layers, which also supports 
the argument that trees at the OS have access to and are 
using water from deeper in the soil profile [Bowling et ai , 
2003]. 

[ 54 ] A comparatively shallower rooting depth parameter 
value for the YS is reasonable given the site’s stage of 
development following disturbance [Law et ai , 2001c]. A 
similarly shallow rooting depth parameter value, though, for 
the MS raises questions regarding whether the root system 
has had adequate time to develop or whether there are site 
factors that could be limiting root system development. 
Ordinarily, individual ponderosa pine trees would be 
expected to establish their root systems within a few years, 
and pure stands would have fully developed root systems 
within a few decades [Oliver and Ryker , 1990; Stone and 
Kalisz , 1991]. An analysis of empirical height-age relation- 
ships similar to that used to develop site index curves for 
forestry applications [e.g., Meyer , 1938] suggests that both 
the MS and YS may be poorer quality sites compared to the 
OS, which might reflect possible differences in rooting 
depths at the sites. However, these site differences are 
probably not nutrient related. Recently, Kelliher et ai 
[2004] reported that concentrations of total carbon (C), total 
nitrogen (N), and mineral N in the top 10 cm of soil, where 
most fine roots are located, were similar at YS and OS. At 


deeper depths, though, the YS had higher concentrations of 
total C and total N but not mineral N. Nitrogen minerali- 
zation rates at the two sites were also similar, and they 
concluded that mineralization rates of both soil and litter C 
at the sites was limited more by the availability of water 
than by N [Kelliher et ai, 2004]. It is possible, though, that 
some physical obstruction, perhaps related to site geology 
or topography, may be restricting root system development 
at the MS. 

[ 55 ] In general, our simulations support the view that 
older forest stands, with well-established root systems, are 
better able to buffer effects of intra-annual and interannual 
variation of climate than younger stands with less-devel- 
oped root systems that grow in semi-arid regions like central 
Oregon. One concern of a potential change in climate is that 
there will be increased variation in individual extreme 
precipitation events as well as increased variation in total 
annual precipitation [Easterling et ai, 2000]. Our simula- 
tions suggest that younger forests, under circumstances of 
increased climatic variability, would be subject to higher 
interannual variation of gross carbon uptake and be more 
susceptible to climatic extremes. In a recent review, Weltzin 
et ai [2003] voiced similar concerns about the potential 
effects of changes in precipitation regimes on ecological 
systems. In addition, much of the western United States 
(including semi-arid regions with ponderosa pine) is sub- 
ject to both natural (e.g., wildfire) and human-caused (e.g., 
logging) disturbances, which have a direct impact upon the 
stage of development and structure of forested ecosystems. 

Thus, in addition to the ecological implications of our 
results for ecosystem functioning, our simulation results 
may have further implications for forest management 
policies. 

3.6. Implications for Net Ecosystem Exchange 

[ 56 ] The model-based annual estimate of NEE for the YS 
in 2001 compared favorably with the estimates derived from 
eddy flux data (Table 5). However, for 2002, modeled NEE 
was considerably larger than the flux-based estimates. At 
the MS, the larger modeled NEE (compared with measured) 
is related to the much higher simulated GPP at the site. At 
the YS, it is unclear why modeled NEE was higher than 
measured NEE in 2002, but it may be related to a higher 
estimate of foliage respiration (R/) in 2002 than 2001, 
although R e is similar between the two years (Table 7). 

[ 57 ] At the OS, NEE was lowest in 2000 and highest in 
1999, but among the 4 years, the site was a consistent and 
stable net carbon sink (Table 7). In contrast, NEE at the YS 
varied widely. In 1999, the site was a small net carbon 
source but was a net sink for carbon thereafter. Similarly, 
the 2 years of available data for the MS suggest that NEE at 
the site was highly variable and that the site was a net sink 
for carbon. Across the 4 years, interannual variation of C0 2 
fluxes from the soil surface ( R s ) was similar among the three 
sites and accounted for >70% of annual R e at the YS and OS 
and ~64% at the MS. Consequently, interannual variation 
of R e at the sites was controlled mostly by interannual 
variation of R s . At the YS, the variation of R e was also 
controlled in part by interannual variation of Rf as leaf 
area increased. Because GPP at the OS was stable over the gg 


14 of 17 


GB4007 


SCHWARZ ET AL.: CARBON AND WATER FLUXES IN PONDEROSA PINE ECOSYSTEMS 


GB4007 


Table 7. Annual Estimates of Net Ecosystem Exchange (NEE) Determined From Simulations of GPP and Annual Estimates of 


Respiration Derived From Measurements at Each of the Study Sites 3 





YS 





MS 





OS 



R, 

R r 

Re 

GPP 

NEE 

Rs 

R ( 

Re 

GPP 

NEE 

R , 

Rf 

Re 

GPP 

NEE 

1999 

453 

111 

590 
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“Units are g C m 2 yr *. R s , C0 2 efflux from the soil surface; foliage maintenance respiration; R e , whole ecosystem respiration, R e = R s + R/ + 
sapwood respiration (R w ) + annual decomposition of woody material (see section 2.8). NEE = GPP - R e . Dashes indicate no data available. 


4 years, interannual variation of NEE at the site was 
determined mostly by the interannual variation of R e , 
whereas the interannual variation of GPP probably had a 
much bigger effect on the interannual variation of NEE at 
the YS and MS (Table 7). 

[58] Determining the sign, magnitude, and interannual 
variation of NEE in various ecosystems is critical to 
developing a better understanding of global carbon budgets 
[Schimel et al. , 2001]. Our results suggest that the controls 
on the interannual variation of NEE in ponderosa pine 
forests, even with similar climate, are not simple and 
consistent across age classes but may differ according to 
development stage. Although the robustness of the NEE 
calculations using our approach is unclear, the apparent 
transition of the YS from a net carbon source to a net carbon 
sink highlights the importance of forest recovery following 
disturbance in determining terrestrial carbon balances [Law 
et al ., 2001c; Schimel et al , 2001; Goodale et al. , 2002]. 

[ 59 ] Our results contribute to the debate regarding the 
relative importance of R e versus GPP in controlling the 
interannual variation of NEE. Some studies of forests 
suggest that R e varies more than GPP on an annual basis 
[Valentini et al. , 2000; Janssens et al , 2001; Valentini et al . , 
2003]. Other studies indicate that GPP varies more than R e 
[Arain et al , 2002; Aubinet et al. , 2002], while other studies 
show that the interannual variation of R e and GPP is roughly 
the same [Barford et al , 2001; Suni et al , 2003; Wang et 
al. , 2004]. Our results are important because they suggest 
that in forests with similar climate, the relative importance 
of R e and GPP may depend upon the development stage of 
the forest. 

4. Conclusions 

[ 60 ] Meteorological data collected at the study sites indi- 
cated that precipitation was the dominant source of climatic 
variation over the 4-year period, 1999-2002, and that the 
variation in annual precipitation exceeded 50%, while 
summer precipitation varied more than six-fold. Model 
simulations at the three sites suggested that interannual 
variation of GPP over the 4 years ranged from 7% (OS) 
to >50% (MS). Simulations that assessed the effects of 
seasonal drought at the sites suggested that the YS and MS 
were relatively more constrained by seasonal drought stress 
than the OS throughout the 4-year period, such that inter- 


annual variation of climate had a larger effect on the YS and 
MS and very little effect on GPP at the OS, which was 
buffered by deeper rooting. At the rapidly developing YS, 
changes in stand structure, such as increasing leaf area 
associated with vigorous growth, appear to have larger 
effects on carbon and water fluxes than variation in climate, 
although effects of these changes may interact with other 
biotic effects including shifts in carbon allocation and 
whole-tree hydraulic conductance. Interannual variation of 
NEE, based on simulations of GPP and empirical estimates 
of R e , was also less variable at the OS than at the other two 
sites, and whereas interannual variation of NEE at the OS 
during the 4 years was probably controlled mostly by the 
variation of R e , interannual variation in NEE at the YS and 
MS appeared to be more strongly controlled by the inter- 
annual variation of GPP. 

[6i] Our results suggest that the interannual variation of 
precipitation is probably the dominant control on carbon 
and water vapor fluxes in temperate coniferous forests 
growing in semi-arid regions. Additionally, older forest 
stands with well-established root systems appear to be better 
able to buffer the effects of both seasonal drought stress and 
interannual climatic variation than younger stands. In many 
forested landscapes, mosaics of different forest development 
stages are present because of the effects of wildfire and 
timber harvesting (Figure 1). Carbon uptake and water 
vapor exchange among stands with different structural 
characteristics and stages of development can show varying 
responses to the interannual variation of climate, and the 
aggregate pattern of carbon and water fluxes at the land- 
scape level will likely depend upon the distribution of 
different development stages. If there is a higher proportion 
of younger stands in the landscape because of disturbances 
like intensive harvesting or wildfire, then GPP across the 
landscape is likely to undergo large interannual variations in 
response to climate, whereas a higher proportion of old 
stands would probably dampen the interannual fluctuation 
of GPP. In western North America where ponderosa pine 
forests are common, there are relatively few old stands; thus 
the interannual variation of GPP in this region is likely to be 
large. Because interannual variation of NEE in younger pine 
forests is probably more strongly related to the interannual 
variation of GPP than R e , NEE across the region is also 
likely to vary considerably in response to the interannual 
variation of climate. 
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[i] This study evaluates the prediction of heat and moisture fluxes from a new land 
surface scheme with eddy correlation data collected at the old aspen site during the Boreal 
Ecosystem- Atmosphere Study (BOREAS) in 1994. The model used in this study couples a 
multilayer vegetation model with a soil model. Inclusion of organic material in the 
upper soil layer is required to adequately simulate exchange between the soil and 
subcanopy air. Comparisons between the model and observations are discussed to reveal 
model misrepresentation of some aspects of the diurnal variation of subcanopy 
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1. Introduction 

[ 2 ] Modeling approaches for surface fluxes of heat, mois- 
ture and carbon dioxide can be classified into single-source 
models, two-source models and multilayer models. In the 
single-source model, evaporation is determined as if the plant 
canopy were a partly wet plane at the lower boundary of the 
atmosphere using bulk aerodynamic resistance and stomatal 
resistance. The bulk stomatal resistance in the single-source 
model is less well behaved than leaf stomatal resistance in 
two-source or multilevel models since it is not a purely 
physiological parameter [Raupach and Finnigan , 1988]. 
Results from single-source models can be particularly sensi- 
tive to the roughness lengths for heat and moisture, which can 
behave erratically over vegetated surfaces. 

[ 3 ] In two-source models [ Kustas , 1990; Sellers et al, 
1986; Choudhury and Monteith, 198 S; Norman etal, 1995], 
an explicit single vegetation layer is considered separately 
from the ground surface. Although this model is more 
realistic than the single-source model, vertical structure of 
the canopy is not resolved. In actual canopies, the stomatal 
resistance depends significantly on height within the canopy 
because radiation, turbulence transfer and water supply from 
the root system vary with height. Consequently, multilayer 
models have been developed to simulate tall canopies such as 
forests [*Sm etal., 1996; Albertson etal , 2001], These models 
aim to describe not only the evaporation from the entire 
canopy, but also the partitioning of the evapotranspiration 
between various parts of the canopy together with other 
aspects of the canopy microclimate such as profiles of leaf 
and air temperature and air humidity. The price of such details 
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is more complexity and parameter input requirements. The 
choice of model complexity depends on the purpose and 
availability of computer resources. 

[4] Even when multilevel canopy models approximate 
the actual canopy processes, the coupling between the 
atmospheric boundary layer and canopy models requires 
estimation of aerodynamic quantities [Ska et al, 1999] and 
compliance with Monin-Obukhov similarity theory. In the 
roughness sublayer immediately above the canopy, the flux- 
gradient relationship based upon Monin-Obukhov similarity 
theory can significantly underestimate scalar fluxes 
[Simpson et al. , 1998; Kaimal and Finnigan , 1994]. Param- 
eterizations for the roughness sublayer [Cellier and Brunet , 
1991; Wenzel et al , 1997; Physick and Garratt , 1995] are 
difficult to verify from observations, partly because of the 
potentially large horizontal gradients on the scale of the 
roughness elements [e.g., Katul et al, 1999]. Such micro- 
scale heterogeneity can contribute to vertical flux diver- 
gence due to increasing footprint with height. Two-source 
models apply Monin-Obukhov similarity theory by assum- 
ing that the aerodynamic temperature for Monin-Obukhov 
similarity is equal to the canopy air temperature. Although 
this simplification can lead to significant errors [Sun et al , 
1999], alternative procedures have not been developed and 
we employ the same approximation in this study. 

[5] The canopy turbulence is driven partly by coherent 
eddies of canopy scale sometimes leading to locally coun- 
tergradient fluxes [Raupach et al , 1996]. Existing numer- 
ical models of canopy turbulence span a wide range of 
complexities from K theory to higher order closure model- 
ing and LES modeling [Katul and Albertson , 1998; Shaw 
and Schumann , 1992]. Although K theory does not describe 
countergradient fluxes, for practical purposes, K theory 
remains an adequate approximation for some applications 33 
[Dolman and Wallace , 1991] and is still used in many large- 
scale models [Bonan, 1996; Cotton et al , 2003]. 
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[6] The coupling between the canopy air and soil through 
fluxes of heat and moisture can be governed by a shallow soil 
layer of high organic content above the mineral soil [Pauwels 
and Wood , 1999a; Van de Wiel etal ., 2002]. The organic layer 
is not usually included in soil models. Blanken et al [1997] 
and others have noted that leaf litter on the forest floor can act 
as a quasi-insulator and therefore promote large subsurface 
temperature gradients in addition to greater daytime warming 
of the subcanopy air in the pre-leaf-out period. More specif- 
ically, the soil surface layer with high organic content is 
normally characterized by smaller thermal conductivity and 
high porosity compared to mineral soils [Letts et al, 2000]. 
Pauwels and Wood [ 1 999a, 1 999b, 2000] also pointed out the 
importance of an organic surface layer and included an 
organic surface layer in TOPLATS. The model performed 
well in simulating surface fluxes above the canopy but 
slightly overestimated the ground heat flux during the course 
of the day [Pauwels and Wood , 1999b] in BOREAS [Sellers 
et al., 1995]. Although the modeled ground heat flux is 
sensitive to parameters such as moss thickness, thermal 
conductivity and heat capacity [Pauwels and Wood , 
1999b], subcanopy turbulence also plays an important role 
in determining ground heat flux. In this study, we focus on 
evaluation of the subcanopy flux and interaction between the 
surface organic layer and the subcanopy air. 

[ 7 ] To examine the interaction between microclimate and 
physiology, we have coupled the multilevel canopy model 
of Williams et al [1996] with a soil model and atmospheric 
boundary layer model for use in regional models. This study 
evaluates the offline performance of the canopy and soil 
models in terms of vertical structure within the canopy and 
interaction with an organic surface layer. The model will be 
compared with data collected at the old aspen site during 
BOREAS. 


2. Model 

[8] The land surface model is based on the canopy model 
(soil plant atmosphere (SPA) model) of Williams et al 
[1996, 2001], coupled to a multilayer soil model with snow 
and frozen soil physics [Mahrt and Pan , 1984; Koren et al, 
1999; Peters-Lidard et al, 1998] and the surface runoff 
scheme of Schaake et al [1996]. The canopy model 
computes the stomatal resistance in each canopy layer to 
maximize daily carbon gain per unit leaf nitrogen content, 
within the limitation of canopy water storage and transport 
of water from soil to the canopy [Williams et al, 1996]. The 
radiation routines model the incidence, interception, absorp- 
tion and reflectance of PAR (Photosynthetic Active Radia- 
tion), near infrared radiation (NIR) and longwave radiation 
in each canopy layer [Amthor, 1994; Amthor et al, 1994]. A 
spherical leaf angle distribution is assumed. 

[ 9 ] Here we apply two vegetation layers with the same 
thickness. Four soil layers are employed with thicknesses of 
0.1, 0.3, 0.6 and 1 m, where the top layer includes organic 
material. We refer to the above papers for a description of 
the model components and here describe only parts of the 
model where changes are made. 

2.1. Subcanopy Processes 

[ 10 ] The total surface fluxes are partitioned into vegeta- 
tion and ground components. Transport is formulated in 
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Figure 1. The transfer diagram for the canopy system. 


terms of the atmospheric eddy diffusivities and leaf boundary 
layer resistances (Figure 1). The subcanopy air temperatures 
(Tj and T 2 ) and the water vapor mixing ratios ( q x and q 2 ) 
are calculated using the leaf temperature (7J V ), the saturation 
mixing ratio ( q iv ) evaluated at the leaf temperature, the ground 
temperature ( T g ), the ground saturation mixing ratio (q g ) 
for the ground temperature and resistances by solving the 
following equations for water vapor and temperature by 
iteration: 
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where q a and T a are the mixing ratio and the air temperature 
at the reference height (39 m from the ground) above the 
canopy, R a is the aerodynamic resistance in the surface 
layer, R k is the subcanopy resistance, R b is the leaf boundary 
layer resistance and R sun and R sha are the stomatal resistance 
for sunlit and shaded leaf areas, respectively, f w is the 3 g 
wet fraction of leaf and (3 is the soil wetness function 
(equation (13)). L sun is the sunlit leaf area index, L s h a the 
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shaded leaf area index and L and S are the leaf area index and 
stem area index, respectively. While transpiration occurs on 
one side of the leaf, sensible heat flux occurs from both sides 
of the leaf due to mechanical mixing. Therefore a factor of 2 is 
used in equations (3)-(4). The value of 0.93 in equations (3) 
and (4) accounts for the differences in molecular diffusivities 
between heat and water vapor. It is assumed that the air 
adjacent to the shaded leaf surface is the same as that 
adjacent to the sunlit leaf surface due to sufficient mixing. 
Therefore no differentiation of leaf boundary layer resistance 
was made between shaded and sunlit leaf surfaces. The leaf 
temperature is calculated from the leaf energy balance 
equation. Neglecting any metabolic and physical storage, the 
leaf energy balance equation becomes 

R n = H v + \E V ( 1 -f w ) -l- \E c f w , (5) 

where H v is the leaf sensible heat flux, R n is the net radiation 
of the leaf surface, E v the transpiration and E c the direct 
evaporation from wet leaf surfaces. The ground temperature 
is calculated from the ground surface energy balance. 

[i i ] The aerodynamic resistances for momentum and heat 
transfer are calculated using the surface layer similarity theoiy 
for the eddy diffusivity. The subcanopy resistance for the heat 
and moisture fluxes between the soil surface and the lowest 
canopy layer and between canopy layers is estimated as 

Rk = !~K h ■ (6) 

The eddy diffusivity within the canopy layer for temperature 
and moisture is assumed to decrease exponentially from the 
canopy top downward toward the ground surface [ Bonan , 
1996] 

Kh — Ksfc exp(-a(l - ^))i (?) 

where K sfc is the eddy diffusivity at the canopy top 
calculated from Monin-Obukhov similarity, a is a non- 
dimensional constant, z is the height above ground and H is 
the canopy height. 

[12] The leaf boundary layer resistance is calculated as in 
the work of Jones [1992] 

R b = 100 (d/u) 0! , (8) 

where d is the characteristic dimension of leaf and u is the 
wind speed in each canopy layer. The default value of 
0.08 m from Williams et al [1996] is used for d in this 
study. The wind profile within the canopy is assumed to 
decrease exponentially downward as 

U(z) = t/ c exp(-a(l ( 9 ) 

where U c is the wind speed at the canopy top. 

2.2. Surface Organic Layer 

[13] The topsoil layer is generalized to include organic 
content. The thermal conductivity is computed as 

( 10 ) 


where K t is the thermal conductivity of each material 
component (Table 1) and f is the volume fraction of each 
material. The volumetric heat capacity of the topsoil layer, 
C (J K _1 m -3 ) is represented as 

C = VCifi, (11) 

where C, is the volumetric heat capacity of the i th soil 
component. For the thermal conductivity of the mineral soil, 
Johansen’s parameterization [Peters-Lidard et al., 1998] is 
used. 

[m] For soil evaporation, we adopted the commonly used 
expression 

E = PE P , ( 12 ) 

where E p is potential evaporation, which is calculated by a 
Penman-based energy balance approach and (3 is calculated 
as 


where © s is saturation soil moisture, 0! is soil moisture at 
first soil layer and 0 W is the wilting point. 

[15] The organic soil has a small bulk density, typically 
about 0.13 g cm -3 , while mineral soil has a typical value of 
about 1.3 g cm -3 for the BOREAS sites [Halliwell and 
Apps, 1997]. The organic material usually includes highly 
permeable fibric peat near the surface [Letts et al ., 2000], 
which dries out quickly, corresponding to rapidly decreas- 
ing hydraulic conductivity and decreasing surface evapora- 
tion. To include this effect into the above formulation, we 
used the following values as the saturation point and wilting 
point. 

e s = 0.87 (14) 

S w = 0.22 (15) 

The saturation point is based on bulk density and wilting 
point is used as a tuning parameter rather than pure physical 
quantity. A similar approach was employed by Pauwels and 
Wood [1999a, 1999b] where the soil resistance was 
calibrated for each tower site to represent the reduction of 
evaporation from a surface moss layer. 

3. Data 

[16] We compared the land surface model with eddy 
correlation data collected at 39 m on the old aspen tower 
in BOREAS 1994 and at 4 m on a small tower approx- 
imately 40 m from the main tower. The study site 
(56.629°N 1 06.200° W) was located in Prince Albert 
National Park, approximately 50 km NNW of Prince 
Albert, Saskatchewan, Canada. The site lies near the 
southern limit of the boreal forest with the transition to 
parkland occurring approximately 15 km to the south- 
west. The soil texture is sandy loam covered by about 4Q 
8 cm of organic material. A natural fire occurred approx- 
imately 70 years ago resulting in an even aged stand of 
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Table 1. Thermal Properties of Soil Constituents From Farouki 
[1986] 


Material 

Heat Capacity, 

xlOV/r'm - 3 

Thermal Conductivity, 
Wm~ x K~ { 

Quartz 

1942 

8.4 

Soil minerals 

1942 

2.9 

Soil organics 

2503 

0.25 

Water 

4186 

0.6 

Ice 

1883 

2.5 

Air 

1.2 

0.026 


aspen with a mean canopy height of 21.5 m and mean 
diameter at the 1.3 m height of 20 cm. Crown space was 
limited to the upper 5-6 m, beneath which was a 
branchless trunk space. The understory was dominated 
by a uniform cover of hazelnut with a mean height of 
2 m. The fetch was at least 3 km in all directions. 

[n] Two periods were selected for model evaluation 
and sensitivity tests during which skies were mostly clear 
except for a few short rain events. The characteristics of 
the two periods are described in Table 2. The first period 
is pre-leaf-out and snow free. According to Blanken et al 
[1997], the leaf-out began in the third week of May. The 
second period is well after full leaf-out when the leaf area 
index (LAI) was approximately constant with time. The 
parameters used in SPA are listed in Table 3. To force the 
model in an offline mode, we used meteorology data 
observed at 39 m on the main tower [Hartog and 
Neumann , 2000], precipitation and longwave radiation 
from the Airborne Fluxes and Meteorology data set 
[Osborne et al, 1998] and observed soil data from Black 
[2000]. The model time step is 15 min. Meteorological 
forcing data are available at 30-min intervals, and pre- 
cipitation and longwave radiation are available at 15-min 
intervals. We have interpolated meteorological forcing 
data into 15-min intervals. 

[is] While we consider the data to be the “truth” for 
evaluation of the model, we must recognize a variety of 
errors. Random flux errors for individual 30-min. records 
may be greater than 10% and often greater than 20% for 
nonstationaiy transition periods and some nocturnal periods 
[Mahrt, 1998]. Subcanopy measurements may be influenced 
by subcanopy heterogeneity. Measurements of wind and 
temperature in the subcanopy are not representative of the 
entire subcanopy layer particularly at night when a strong 
surface inversion generally forms in the lowest 5 m [Mahrt et 
al., 2000]. 

[ 19 ] Nakamura and Mahrt [2001] found that the above- 
canopy eddy correlation measurements in BOREAS were 
generally within the roughness sublayer, below the sur- 
face layer where Monin-Obukhov similarity theory 
applies. Model uses Monin-Obukhov similarity theory, 
which generally underestimates mixing in the roughness 
sublayer. As a probable consequence, the observed verti- 
cal temperature difference between the two observational 
levels (13 m, 39 m) is significantly smaller than 
that predicted using Monin-Obukhov similarity theory 
(Figure 2), especially in stable conditions, implying 
greater mixing compared to Monin-Obukhov similarity 
theory. The influence of the height dependence of the 
footprint on the observed vertical temperature difference 


and fluxes cannot be assessed here, which is always a 
concern in the roughness sublayer. 

4. Comparison With Aspen Data 

[ 20 ] The model was run for two 10-day periods, one in 
the pre-leaf-out and the other in the post-leaf-out period. 

4.1. Pre-Leaf-Out Period 

[ 21 ] Comparison of the 10-day averaged diurnal cycle 
between the model and the observations (Figure 3) indicates 
that the above-canopy and understory fluxes are simulated 
by the model reasonably well during die pre-leaf-out period, 
with some exceptions. Table 4 shows the statistical compar- 
isons between model-derived and tower-observed fluxes. 
Because of the small LAI during this period, the total latent 
heat flux is dominated by soil evaporation. The relatively 
low correlation between the model and observed evapora- 
tion (Table 4) is due to poor model simulation of the direct 
evaporation from the soil immediately after rainfall. 
Although the rainfall amount is small, the observed subcanopy 
moisture flux, including direct evaporation from the soil and 
canopy, reaches about 200 Wm“ 2 immediately after rainfall. 

[ 22 ] The observed sensible heat flux of nearly 1 50 Wm' 2 in 
the subcanopy for the pre-leaf-out period implies significant 
buoyancy-generation of turbulence energy within the subcan- 
opy. In sparse canopies with low LAI, the use of the usual 
constant extinction coefficient (equation (7)) can underesti- 
mate turbulent mixing within the canopy, resulting in larger 
vertical temperature gradients in the subcanopy and higher 
ground surface temperature than observed during daytime. 
The model shows a warm bias in 4-m subcanopy temperature 
during daytime (Figure 3), consistent with underestimation of 
mixing, although advection and errors in the model radiative 
transfer could also be factors. The apparent cold bias in the 
nocturnal subcanopy, could be partly due to the location of the 
observed temperature near the top of a strong surface inversion 
of about 5-m depth. Since the model does not consider 
subcanopy stability in the heat transfer in the subcanopy, the 
model does not capture the strong nocturnal surface inversion, 
which results in overestimation of the understoiy sensible heat 
flux. 

4.2. Post-Leaf-Out Period 

[ 23 ] The post-leaf-out period is characterized by a signif- 
icantly higher LAI and lower soil moisture content than in 
the pre-leaf-out period (Table 2). The simulated latent heat 
flux and soil heat flux are generally close to the observed 
values (Figure 4), although the understory latent heat flux is 
overestimated by 30 Wm -2 during the late afternoon from 
1500 LST to 1800 LST. The transpiration is controlled by 
the leaf boundary layer resistance and the stomatal resist- 
ance. In unstable conditions with significant wind, the leaf 
boundary layer resistance is usually smaller than the sto- 


Table 2. Characteristics of Selected Periods 



PAI(LAI) a 

Initial Topsoil 

Period 

Hazelnut/Aspen 

Moisture 

Pre-leaf-out 5-14 May 

0.3(0.1)/1.2(0.3) 

0.23 

Post-leaf-out 30 July to 8 Aug. 

3.4(3.2)73.2(2.3) 

0.11 


*PAI is the plant area index {the sum of leaf area index and stem area 
index) and LAI is the leaf area index. 
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Table 3. Model Parameters 


Parameter 

Value 

Unit 

Source 

Stem specific hydraulic Conductivity 

20 

m, mol, m -1 , s~\ MPa -1 

estimated 

Root resistivity 

50 

MPa s g mmof 1 

estimated 

Total fine root biomass 

657 

g m -2 

Steele et al [1997] 

Canopy layer capacitance 

8000 

mmol MPa -1 

Williams et al [1996] 

Minimum leaf water potential 

-2.3 

MPa 

Kimball et al. [1997] 

RuBP carboxylation catalytic Rate coefficient at 30°CK c 

37.8 

\i mol g -1 N s' 1 

Williams et al [2001] 

Electron transport rate Coefficient at 30°CK y 

49.0 

\i mol g -1 N s' 1 

Williams et al [2001] 

Soil color 

4 

no unit 

Kimball et al [1997] 


matal resistance. Therefore the transpiration is limited by 
stomatal resistance and is less sensitive to leaf boundary 
layer resistance. However, when the wind speed is weak, 
transpiration can be more limited by the leaf-boundary layer 
resistance. In late afternoon when the subcanopy air 
becomes stably stratified, the wind speed within the canopy 
becomes weak due to less downward transport of momen- 
tum. The overestimation of subcanopy wind speed causes 
underestimation of leaf boundary layer resistance, which 
causes overestimation of transpiration in late afternoon. 

[ 24 ] Model errors for the sensible heat flux are greater 
than those for the latent heat flux. The model overestimates 
the subcanopy sensible heat flux throughout the diurnal 
period (Figure 4 and Table 4), possibly due to failure to 
resolve the vertical structure of the understory. The model 
employs only two layers of vegetation for application in 
regional models. The two vegetation layers, here each 10 m 
deep, assume uniform plant distribution within each layer 
whereas the actual vertical distribution includes a thick 
overstory of 5-6 m thickness and a 2-m dense understory 
near ground. The layer-averaged resistance in the lower 
canopy layer underestimates the leaf boundary layer resis- 
tance for the 2-m understory and therefore overestimates the 
subcanopy sensible heat flux from the understory. Note that 
transpiration is less sensitive to leaf boundary layer resist- 
ance since it is limited by stomatal conductance in unstable 
conditions, compared to sensible heat flux that depends 
more on leaf boundary layer resistance. 

[ 25 ] In the transition periods with low sun angle, the 
modeled subcanopy is stably stratified while the atmosphere 
above the canopy is still unstably stratified. The subcanopy 
mixing in the model depends on the stability above the 
canopy and therefore does not recognize the difference 
between the stability above and within the canopy. As a 
result, the model incorrectly predicts large negative sensible 
heat flux in the subcanopy during the transition periods 
whereas the observed subcanopy heat flux is small (Figure 4). 

[ 26 ] The model overestimates the daytime sensible heat 
flux above the canopy by 50 W m -2 . This overestimation 
could be due to omission of the heat storage within the 
canopy and advection of temperature. 

[ 27 ] The simulated soil heat flux agrees with the obser- 
vations reasonably well during the day, but is overpredicted 
at night (Figure 4). The modeled temperature at the 4-m 
level is in good agreement with the observed temperature 
during the day but it is too cool at night giving larger 
temperature gradient than that observed in the upper canopy. 
The apparent cold bias could be due to location of the 
observed temperature just above the subcanopy surface 
inversion (section 4.1) or could be due to model overesti- 
mation of the radiative cooling of leaves. 


[ 28 ] The leaf radiative cooling is coupled to the leaf 
boundary layer resistance (equation (8)) through its influ- 
ence on leaf temperature in the leaf energy balance. The leaf 
boundary layer resistance represents the transfer of heat and 
moisture between the leaf surface and adjacent air and in the 
model is only a function of wind speed. In unstable 
conditions, the mixing in the canopy is significant, such 
that the canopy wind speed and air temperature are close to 
those adjacent to the leaf surface. However, in stable 
conditions, large differences exist between average canopy 
wind speed and air temperature and wind speed and air 
temperature adjacent to the leaf surface due to suppression 
of mixing by the stable stratification. Without consideration 
of these differences, the heat exchange between the leaf and 
the atmosphere is overestimated in stable conditions, which 
reduces the temperature difference between the radiatively 
cooled leaf surface and the air. The resulting overprediction 
of leaf temperature causes more emitted longwave radiation 
and therefore larger net radiative cooling of the canopy. Due 
to the overcooling in the subcanopy, the subcanopy air 
temperature decreases to the saturation point causing over- 
estimation of understory condensation (Figure 4). 

5. Influence of Organic Material 

[ 29 ] To examine the effect of the surface organic layer on 
the subcanopy processes, simulations were performed with 
and without organic material in the upper soil layer. During 



z/L 


Figure 2. The modeled vertical temperature difference 
(r I3m — T 39m ) minus the observed vertical temperature 42 
difference between the displacement height (13 m) and 
reference height (39 m) as a function of z/L. 
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Sensible heat flux (understory) 


Latent heat flux above canopy 





Figure 3. Composite diurnal variation for the observed (circles) and modeled (solid line) variables for 
the pre-leaf-out period. 


the pre-leaf-out period, the simulation with no organic 
material substantially overestimates soil evaporation and 
heat flux into the ground and underestimates understory 
sensible heat flux (Figure 5). The simulation with organic 
material predicts fluxes much closer to the observed values. 
The unrealistically large soil evaporation in the simulation 
with no organic material leads to rapid drying of the soil. 
The organic material acts as a partial insulator between the 
soil surface and deep soil. As a result, the surface organic 
layer makes the subcanopy more unstable during daytime 
and more stable during night. This improves the comparison 
with the observations. 

[ 30 ] The impact of the organic material is smaller after 
leaf-out because the understory reduces exchange between 


the ground and the atmosphere (Figure 6). The overestima- 
tion of soil evaporation in the simulation without an organic 
layer leads to underestimation of the vapor pressure deficit 
in the subcanopy, which in turn leads to underestimation of 
transpiration. Therefore the sensitivity of latent heat flux to 
the organic layer is reduced. The simulation with no organic 
material also overestimates the downward ground heat flux 
during the daytime, which results in smaller sensible heat 
flux in the subcanopy. 

6. Conclusions 

[ 31 ] The land surface scheme tested in this study roughly 
captures the main energy partition between the understory 


Table 4. Results of Model Application® 


Period 

Variable 

Mean Observation, W m 2 

Mean Simulation, W m 2 

R, Dimensionless 

RMSE, W m' 2 

Pre-leaf-out 

H u 

40.62 

40.65 

0.96 

33.62 


H a 

96.51 

71.76 

0.95 

54.71 


LE U 

11.37 

12.69 

0.62 

15.48 


LE a 

17.06 

19.76 

0.59 

21.62 


G 

-12.36 

- 9.00 

0.92 

14.35 

Post-leaf-out 

H u 

2.39 

-0.51 

0.69 

11.48 


H a 

10.15 

4.63 

0.88 

35.30 


LE U 

21.46 

23.76 

0.81 

23.6 


LE a 

78.02 

74.49 

0.93 

36.42 


G 

-5.03 

-2.52 

0.95 

5.29 


"Subscripts u and a represent understory and above the canopy respectively, R is the correlation coefficient and RMSE is the root mean square error. 
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Latent heat flux (understory) Sensible heat flux (understory) 







Figure 4. Composite diurnal variation for the observed (circles) and modeled (solid line) variables for 
the post-leaf-out period. 


Latent heat flux (understory) Ground heat flux 





Figure 5. Comparison between the model with (solid line) and without (dash-dot line) the organic layer 
for the pre-leaf-out period and the observed values (circles). 
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Figure 6. Comparison between the model with (solid) and without (dash-dot line) the organic layer for 
the post-leaf-out period and the observed values (circles). 


and overstory during both pre-leaf-out and post-leaf-out 
periods. The model also partitions the net energy into latent 
heat flux and sensible heat flux reasonably well, provided 
that the modeled upper soil layer includes organic material. 
Failure to include such material causes substantial overes- 
timation of soil evaporation and subsequent errors in the 
subcanopy fluxes, particularly prior to leaf-out. The organic 
material reduces the heat exchange between the ground 
surface and deeper soil, which makes the subcanopy more 
stable at night and more unstable during the daytime. 

[ 32 ] Because of complex interactions between different 
components of the canopy-ground system, isolating the 
cause of other errors is difficult without detailed observa- 
tions of individual processes. For example, underestimation 
of the leaf boundary layer resistance at night causes over- 
estimation of the transfer of heat from the canopy air to the 
leaf surface, which leads to overestimation of leaf temper- 
ature and overestimation of radiative energy loss from the 
canopy system. Other errors include overestimation of 
mixing in the understory due to failure to resolve the 
vertical structure of the understory. 

[ 33 ] The model does not directly include the important 
influence of diurnal variation of subcanopy stability on the 
mixing. This omission appears to cause overestimation of 
exchange between the ground and subcanopy air in transi- 
tion periods and at night when an inversion layer forms in 
the subcanopy. We are currently examining this influence 
using data from a variety of open canopies. 
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Abstract 


In open canopies, the within-canopy flux from the ground surface and 
understory can account for a significant fraction of the total flux above the 
canopy. This study incorporates the important influence of within-canopy 
stability on turbulent mixing and subcanopy fluxes into a first-order clo- 
sure scheme. Toward this goal, we analyze within-canopy eddy-correlation 
data from the old aspen site in the Boreal Ecosystem - Atmosphere Study 
(BOREAS) and a mature ponderosa pine site in Central Oregon, USA. A 
formulation of within-canopy transport is framed in terms of a stability- 
dependent mixing length, which approaches Monin-Obukhov similarity the- 
ory above the canopy roughness sublayer. 

The new simple formulation is an improvement upon the usual neglect of 
the influence of within-canopy stability in simple models. However, frequent 
well-defined cold air drainage within the pine subcanopy inversion reduces 
the utility of simple models for nocturnal transport. Other shortcomings of 
the formulation are discussed. 


Key words: Canopy mixing, Cold air drainage, Similarity theory, Sub- 
canopy turbulence 

1 Introduction 

The need to understand and quantify biosphere-atmosphere exchange has led 
to substantial interest in within-canopy turbulence. Coherent eddies encom- 
passing the entire canopy often dominate the vertical transport and can lead 
to fluxes counter to the local vertical gradient within the canopy (Raupach 
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et al., 1996, Blanken et al., 1998, Finnigan, 2000). The primary source of 
canopy turbulence is often thought to be downward transport of turbulence 
from above the canopy. 

Modelling within canopy transport normally assumes a Lagrangian or 
Eulerian approach. This study follows the Eulerian approach. A common 
modelling approach extrapolates the mixing coefficient from the surface layer 
above the canopy downward through the canopy using a specified functional 
dependence on height (e.g., Shuttleworth and Wallace, 1985; Bonan, 1996). 
Modelling the flux at a given level includes first-order closure (e.g., Wilson et 
al., 1998), K-epsilon (Katul et al., 2004), second-order closure (e.g., Poggi et 
al., 2004) and large-eddy models (e.g., Shaw and Shumann, 1992, Albertson 
et al., 2001). These models generally require a flux-gradient approximation 
at some level of parameterization, often expressed in terms of a mixing length 
(e.g., Raupach et al., 1996; Wilson et al., 1998; Katul et al., 2004; Poggi et 
al., 2004). 

Using a mixing-layer analogy, Raupach et al. (1996) introduced the shear- 
length for describing turbulence above the canopy top. This length scale has 
been incorporated into several mixing-length formulations. For example, Wil- 
son et al. (1998) proposed a height-dependent mixing length that approaches 
the shear length scale near canopy top. Poggi et al. (2004) employed a length 
scale related to element diameter and the distance from the ground in addi- 
tion to the shear-length scale in the mixing length formulation. On the other 
hand, Massman and Weil (1999) focused on the influence of local leaf area 
index on the mixing length and related the mixing length inversely to the 
local drag coefficient such that the mixing length was smaller at levels with 
more local leaf area index. 

Less attention has been given to the mixing length for momentum and 
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heat within open canopies and the potential role of surface heating and cool- 
ing at the ground/understory surface. Buoyancy effects in the subcanopy 
have been included in a few LES models (e.g. Dwyer et al., 1997; Albertson 
et al., 2001). Brunet and Irvine (2000) investigated the influence of sta- 
bility on the shear-length scale for an open canopy and suggested that the 
primary effect of atmospheric stability is to modify the shear-length scale 
through change of wind speed and wind shear at the canopy top. However, 
they did not examine the influence of local stability on mixing within the 
canopy. Based on eddy-correlation data, Mahrt et al. (2000) found that the 
drag coefficient in the subcanopy at the old aspen site in BOREAS decreased 
with increasing subcanopy stability but did not increase from near-neutral 
to unstable conditions. However, they did not provide a formulation for such 
transport. 

Much of the work above is influenced by observations of the important 
contribution of coherent structures to the within-canopy flux, appearing to 
originate above the canopy, often conceptualized in terms of within-canopy 
sweeps and ejection (e.g. Paw U et al. 1992). Such structures do not pen- 
etrate the strongly stratified subcanopy associated with a open overstrory 
and clear nocturnal conditions (Mahrt, et al . 2000). In open canopies, day- 
time subcanopy turbulence may be generated more by local buoyancy than 
downward transport of turbulence from above the canopy. Resulting under- 
estimation of subcanopy turbulence can significantly influence estimates of 
land-surface exchange in open canopies. 

In this study, we analyze eddy-correlation data from two open canopies 
and construct a simple subcanopy formulation of the mixing length, which 
includes the influence of subcanopy stability. 
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2 Data 


This study analyses eddy-correlation data collected from a mature ponderosa 
pine site in central Oregon, USA (Schwartz et al., 2004). Turbulent fluxes 
were measured above the canopy and at two levels within the canopy, one 
in the crown space (10 m) and one in the trunk space (3 m). The flux at 
3m was measured at a subcanopy tower separated from main tower by 40 
m to avoid disturbance of the flow at the tower base. Wind speed and air 
temperature were measured at 5 levels (3 m, 6 m, 10 m, 20 m, 30 m) on the 
main tower. Air temperature was also measured on the subcanopy tower at 
1 m, 2 m and 3 m. 

This study also analyzes tower data collected at five levels on the main 
tower in the old aspen site in BOREAS (Sellers et al., 1995; Blanken et 
al., 1997). Momentum fluxes were computed from 30-min records. The 
instrumentation at the old aspen site is detailed in Blanken et al. (1997). 
Here we include only the post leafout period when the data were the most 
complete. 

Both canopies are characterized by a LAI of about 3 (Table 1), which we 
refer to as open canopies. While such values are not low compared to sparser 
canopies with a LAI of unity or less. Here, “open” refers to large diurnal 
variation of atmospheric stability within the canopy layer due to radiative 
heating and cooling of the ground/understory surface. The nocturnal tem- 
perature profiles at the old aspen site correspond to a strong surface inversion 
in the lowest 5 m and weaker stratification above (Mahrt et al., 2000). The 
stratification of the daytime aspen canopy layer is near neutral. The diurnal 
variation of the stratification at the pine site is much greater than that at 
the aspen site in spite of similar LAI values for the overstory. The strong 
diurnal variation of stability at the pine site is due partly to more clumping, 
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dryer surface conditions, less understory, higher sun angles, much less cloud 
cover and minimal integrable atmospheric water vapour content, all acting to 
increase the radiative heating and cooling of the ground surface/ understory. 

Fluxes are computed using deviations from 10-minute averages and then 
are averaged over one hour to reduce random flux errors. Fluxes in the sub- 
canopy are subject to some uncertainty from possible flux loss due to path- 
length averaging and expected horizontal heterogeneity above the nonuniform 
understory. Based on our assessment of momentum flux loss due to path- 
length averaging, by comparing with three-dimensional hotfilm anemometry 
(unpublished), we have determined that pathlength average can lead to sig- 
nificant momentum flux loss in the lowest meter above short grass and that 
existing correction formulas are not suitable. However, we have not carried 
out suitable measurements in the subcanopy. 


Table 1 Site description. 


Sites 

Average canopy 
height (m) 

LAI 

displacement 
height (m) 

roughness 
length (m) 

Old aspen 

20.1 

3.0 

13.4 

2.0 

Pine 

15.5 

3.3 

11 

1.2 


The displacement height for the aspen site (Table 1) is taken from Naka- 
mura and Mahrt (2001) while the displacement height for the pine site was 
estimated by trial and error to produce the best fit to Monin-Obukhov sim- 
ilarity theory at the top of the tower for near-neutral conditions. 
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3 Transport in the canopy layer 

3.1 Mixing length 

We express the within-canopy flux in terms of an eddy diffusivity 

W = -K m ^, ( 1 ) 

In the data analysis below, the momentum flux will be computed as the 
magnitude of the vector momentum flux and the mean shear will be computed 
as the magnitude of the vector shear without concern for alignment of the 
wind and shear vectors. 

The eddy diffusivity is formulated as 


— lm.U* 


( 2 ) 


where l m is the mixing length for momentum and n* is the height dependent 
friction velocity in the canopy defined as the square root of the magnitude 
of the momentum flux. The mixing length for momentum 

u, 


< 3 > 

is evaluated directly from eddy-correlation and wind profile measurements at 
the pine and aspen sites. 

The mixing length for momentum at the aspen site was computed from 
the observed fluxes and vertical gradients for three different stability classes, 
based on the stability at the top of the tower (39- m level): the unstable case 
(—1 < (z — d)/L < —0.1), the near neutral case (—0.1 < (z — d)/L < 0.1) 
and the stable case (0.1 < (z — d)/L < 1). For this study, transition periods 
are eliminated and we analyze only daytime data between 1000 and 1600 
local time and nocturnal data between 2100 and 0500 local time. Vertical 
gradients are computed from simple finite differencing. 
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The resulting mixing length for momentum, averaged over all of the quali- 
fying records, increases slowly with height in the subcanopy and more rapidly 
with height above the canopy (Figure 1). The mixing length is much smaller 
in the subcanopy in stable conditions compared to neutral and unstable con- 
ditions, probably due to weaker downward transport of turbulence energy 
into the canopy layer. The profile of the mixing length for the pine site is 
based on only three levels and is not shown. 

To formulate the height dependence of the mixing length, we recognize 
that height above the ground surface influences the within-canopy mixing 
while the mixing length must approach Monin-Obukhov similarity theory in 
the surface layer above the canopy roughness sublayer. For neutral conditions 
within the canopy, we propose the simple height dependence 

^mn " (3z (4) 

where l mn is the neutral value of the mixing length for momentum, z is 

the height above the ground surface. (3 is an empirical coefficient to-be- 

determined and not necessarily equal to the von Karman constant. The 
nondimensional coefficient (3 is estimated from the observed profiles, yielding 
0.18 at the aspen site and 0.30 at the pine site. Either f3 depends on canopy 
architecture, or, it cannot be estimated from our observations within a factor 
of two. 

In order that the mixing length satisfy Monin-Obukhov similarity theory 
in surface layer above the canopy, (3 at the top of the roughness sublayer, z r , 
should approach 

(3 = 0A(z r - d)/z r (5) 

where d is the displacement height and z r can be estimated following Raupach 
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• ■ ■ unstable 
neutral 
■ ■ stable 


Figure 1: Height dependence of the mixing length for momentum estimated 
from vertical profiles at the old aspen site for the post leafout period. 
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(1994) as z r = d+2(h — d) where h is the canopy height. Assigning d = C h, 

(3 = 0.4[1 - t^q\ (6) 

As an exercise, we specify the displacement height to be 2/3h (C=2/3), in 
which case (3 becomes 0.2, which is closer to value at the aspen site. For 
the open canopies of this study, the value of C may be smaller leading to a 
somewhat larger value of (3. Although there is no reason to expect the profile 
fit to Eq. 4 to agree with Eq. 6, we note that they both predict comparable 
magnitudes. Below, we use values of /3 from the profile fit since it avoids 
estimating the height of the roughness sublayer and displacement height. 

3.2 Stability influence on momentum transfer 

In analogy to Monin-Obukhov similarity theory, the stability-dependent mix- 
ing length can be expressed as 


7 - 

'"171 ~~ , ■ 

0m 

Here, the nondimensional wind shear is defined as 

Pzdu/dz 

4>m = • 

u* 

We do not expect (f> m to closely follow Monin-Obukhov similarity theory 
within the canopy because the stress varies rapidly with height. The flux- 
gradient relationship in the subcanopy will be limited in different ways by the 
distance from the ground surface/understory, the spacing between the trunks 
and the distance from the overstory. However, both the heat flux and the 
friction velocity are important influences so that the local Obukhov length, 
A, is one of the governing length scales. So that the stability approaches 
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Monin-Obukhov similarity theory at the top of the roughness sublayer, we 
formulate the subcanopy stability parameter as 


Z.Zr-d 

A 1 * r ' 


( 9 ) 


where the ratio in parentheses allows continuity of the local canopy stability 
functions to that above the roughness sublayer. 

The nondimensional wind shear (Figures 2-3) varies more slowly with sta- 
bility than predicted by the normal stability functions for Monin-Obukhov 
similarity theory at both the aspen and pine sites. The relationship between 
the nondimensional shear and stability exhibits substantial scatter, which 
could be due partly to the difficulty of measuring fluxes in the subcanopy 
(Section 2) and the omission of additional physical influences. In spite of 
the large scatter, inclusion of the influence of subcanopy stability is needed, 
although more data is required to determine the optimum stability function. 
The disagreement between the observations and the usual Monin-Obukhov 
stability functions is greatest at the lowest level of the aspen site, just above 
the understory, where the stability dependence is almost undetectable. Inclu- 
sion of a displacement height for the thick understory of 2-m height increases 
the rate at which the nondimensional shear changes with stability and thus 
increases the stability influence. However, this increase is not nearly large 
enough to allow reasonable approximation with the usual Monin-Obukhov 
stability functions. However, the agreement with the Monin-Obukhov sta- 
bility functions improves with height (Figure 2). 

The smaller observed nondimensional shear (greater mixing) for stable 
conditions, compared to the prediction by the usual stability functions, is 
presumably due to mechanical generation of mixing by the vegetation ele- 
ments, not included in Monin-Obukhov similarity. However, the nondimen- 
sional shear above the canopy is also observed to be smaller than predicted 
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by linear dependence on stability for strong stability (e.g. Beljaars and Holt- 
slag, 1991; Cheng and Brutsaert, 2005). The larger nondimensional shear 
for unstable conditions, most evident at the pine site, could be due to the 
partial constraint of the convective eddies by the height of the overstory, 
which enters as an additional length scale not included in Monin-Obukhov 
similarity. 


3.3 Stability influence on heat transfer 


The pine canopy is characterized by unstable stratification (Figure 4) and 
significant upward heat flux in the daytime and strong stable stratification 
at night. The diurnal variation of the within-canopy stability is less at the 
aspen site as discussed in Section 2. 

The heat transport is expressed in terms of the eddy diffusivity for heat 

ew = -K h ^ (io) 


The eddy diffusivity for heat can be related to the mixing length for heat. 
The mixing length for heat for near neutral conditions is difficult to estimate 
from data because of large scatter. We equate the neutral values of the 
mixing lengths for heat and momentum but allow them to have different 
stability functions. 

The nondimensional temperature gradient is defined as 


f3zd0/dz 


where 


ea- 


rn'd' 

u„ 


( 11 ) 


( 12 ) 


was evaluated at the two eddy-correlation levels at mature pine site, one 
in the crown space (10m) and the other in the trunk space (3m). This 
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Figure 2: The nondimensional shear as a function of stability for 3 different 
levels at the Aspen site and the nondimensional temperature gradient at 
the 18-m level. The top of the canopy is approximately 20 m. The solid 
line corresponds to the stability function from Dyer (1974) applied to 7 z/L, 
where z is the height above the ground surface. 7 = 


13 






| (a) 10m (Crown space) (b) 3m (Trunk space) 



Figure 3: The nondimensional shear and nondimensional temperature gra- 
dient for two levels at the pine site. The top of the canopy is approximately 
15 m. The solid line corresponds to the stability function from Dyer (1974) 
applied to 7 z/L. 7 = 
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relationship was also evaluated at the 18-m level at the aspen site where the 
heat flux and vertical temperature gradients are best defined for evaluation 
of the relationship. Evaluation of vertical gradients must recognize that the 
time-averaged temperature can vary horizontally in the subcanopy on the 
microscale, especially close to the ground surface at the pine site. We use 
only data where the horizontal temperature difference at 3m between two 
subcanopy towers and the main tower is less than 0.5 K, which eliminates 
some of the most stable cases. 

At the pine site, the nondimensional temperature gradient more closely 
follows the Monin-Obukhov stability function compared to that for momen- 
tum (Figure 3). Apparently, the influence of diabatic heating in the canopy 
layer and associated height-dependence of the heat flux does not strongly 
modify the dependence of the flux-gradient relationship on stability. How- 
ever, the influence of stability on the nondimensional temperature gradient 
is not as well defined for the aspen site (Figure 2). At both sites, the scatter 
is large for stable conditions, indicating large errors in the flux calculations 
and/or additional influences on the nondimensional temperature gradient, 
not included in z/A. 

Although the nondimensional temperature gradient for heat exhibits a 
stronger dependence on within-canopy stability compared to the nondimen- 
sional shear, the present data do not allow confident assessment of the ver- 
tical structure of the differences between the mixing lengths for heat and 
momentum. 

With the above findings and uncertainties, we formulate the eddy diffu- 
sivity for heat as 

K h = (13) 

0/1 

where l mn is the mixing length for momentum for neutral conditions (Eq. ??), 
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m* is the height-dependent friction velocity and <ph is the nondimensional tem- 
perature gradient with the usual Monin-Obuhkov stability functions, which 
formally enter here as the stability dependence of the mixing length. Since 
the dependence of the nondimensional shear on stability is less certain than 
that for temperature, we choose to follow the traditional top down approach 
of extrapolating the stress downward from the surface layer using a profile 
model. Ultimately, a separate equation for the momentum budget is needed 
for the canopy layer. 

4 Influence of LAI on canopy momentum flux 

Estimation of the diffusivity for heat from the stability-dependent mixing 
length in the previous section requires information on the mean profile of the 
local friction velocity in the canopy, which we now parameterize based on 
the observed behavior. The ratio of the subcanopy friction velocity to that 
above the canopy increases with increasing instability (Figure 5), presumably 
due to greater mixing in the daytime as suggested by larger within-canopy 
mixing lengths in the daytime (Figure 1). We will include the influence of 
the vertical structure of the LAI on the height dependence of the friction 
velocity. For the aspen site, the downward accumulated Leaf Area Index 
( LAI acc ) is estimated based on the canopy structure of the aspen (Blanken 
et al. 1997). We then relate the momentum flux in the mean wind direction 
to stability and canopy architecture as 

ul = ul 0 exp[-2a{LAI acc f h ] (14) 

where a is a stability-dependent coefficient evaluated below and u* 0 is the 
friction velocity in the surface layer above the canopy. The exponential form 
allows the momentum flux to approach the surface layer value above the 
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canopy. The square root dependence on the accumulated LAI was determined 
as a reasonable fit to near-neutral conditions. Although Eq. 14 adequately 
simulates the vertical profile of local friction velocity at the aspen site, the 
significantly deviations of the observed from the predicted values at the pine 
site (Table 2) either indicate inadequate vertical resolution of the flux profiles 
at the pine site or imply that the local friction velocity for near-neutral 
conditions are affected by factors other than accumulated LAI. The pine site 
is characterized by greater foliar clumping and clumping of trees, allowing 
easier vertical exchange and greater within-canopy momentum flux compared 
to Eq. 14. 


Table 2. Vertical profile of the height-dependent friction velocity nor- 
malized to the value above the canopy for neutral conditions. The ratio is 
defined as the local u * divided by friction velocity above the canopy, u* 0 . 


Aspen site 

Pine site 

z (m) 

LAI acc 

Ratio 

z(m) 

LAI acc 

Ratio 

observed 

estimated 

observed 

estimated 

28 

0 

i 

i 

30 

0 

i 

i 

18.6 

1. 

0.49 

0.37 

10 

1.5 

0.62 

0.29 

5.8 

3. 

0.19 

0.18 

3 

3 

0.13 

0.18 

0.45 

6.5 

0.05 

0.08 






Using the observed fluxes, we have computed a from Eq. 14. We have 
visually fit the stability dependence of a for unstable conditions as (Figure 

5 ) 


a (1-42/L) 0 - 25 ' ^ 

Application of this stability dependence in Eq. 14 leads to reasonable pre- 
diction of the ratio of the subcanopy stress to the stress above the canopy for 
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unstable conditions at both the aspen and pine sites. For stable conditions, 
a appears to be independent of stability (Figure 5). While the above inclu- 
sion of differences between stable and unstable conditions is rather crude, it 
is expected to be an improvement upon the usual practice of neglecting the 
influence of stability on the subcanopy transport. However, with the frequent 
occurrence of drainage flow and capping inversion at the pine site (Section 
5), the vertical communication is too weak to apply Eq. 14 and a separate 
time-dependent momentum equation is needed near the surface. 

5 Nocturnal structure of the pine subcanopy 

Part of the large scatter for the nocturnal flux-gradient relationship at the 
pine site is due to cold air drainage and wind directional shear. While some 
cold air drainage also occurred at the aspen site (Lee, 1998), it seemed to 
be less common and have less impact on the vertical structure of the flow 
within the canopy compared to the pine site. With cold air drainage at the 
pine site (westerly flow near the surface), the nocturnal stratification is much 
stronger and concentrated at the top of the cold air drainage, between 10 and 
20 m (Figure 4). For individual hours, the stratification is often confined to 
a single 5-m observational layer but is vertically spread in the composited 
profile of Figure 4. As a result, the nocturnal vertical temperature gradients 
at the 10 m level are probably often seriously underestimated. A network 
of wind and temperature measurements indicates cold air drainage down the 
eastward descending slope, just west of the tower site. In the absence of 
cold air drainage, the subcanopy flow is generally from the south and the 
stratification above the under story is weaker. 

The joint frequency distribution between the wind speed at 30 m and the 
temperature difference in canopy between 3m and 20m (Figure 6) indicates 

19 


66 


39m 


(a) Pine 10 m (b) Pine 3 m 






Figure 5: The ratio of friction velocities as a function of stability where L a j c 
is the Obukhov length based on fluxes at the 39-m level, assumed to be in 
the surface layer. The solid line represents the prediction of the ratio based 
on the fit to Eq. 15 for unstable conditions and independent of the stability 
for stable conditions. 
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Figure 6: Joint frequency distribution of the 30 m wind speed and the canopy 
inversion strength between 3 and 20 m (DT) at the pine site. Color pattern 
is only for visualization. 
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Figure 7: Joint frequency distribution of subcanopy 3-m wind direction and 
the canopy inversion strength between 3 and 20 m (DT) at the pine site. 
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that most of the nights with inversions greater than a few degrees occur with 
weak winds above the canopy and that strong winds generally correspond to 
weaker inversions of less than a few degrees. Almost all of the strong inver- 
sions occur with westerly subcanopy drainage flow. Weak canopy inversions 
correspond to southerly subcanopy flow more aligned with the flow above the 
canopy (Figure 7). The systematic development of subcanopy drainage flows 
and associated directional shear and strong inversions at the pine site prob- 
ably account for much of the large scatter in the flux-gradient relationship. 
Inclusion of the directional shear associated with the drainage flow requires 
a separate momentum equation for the subcanopy. 

6 Summary 

The structure of turbulent transport within the canopy has been examined 
with eddy-correlation measurements from a relatively open ponderosa pine 
site and from an aspen site. The dependence of the nondimensional temper- 
ature gradient and mixing length for heat on stability approximately follow 
MO similarity theory even though the assumptions required for such simi- 
larity theory are not met. The nondimensional shear and mixing length for 
momentum deviate more from MO similarity. Deviations are most significant 
close to the top of the understory where the mixing coefficient depends much 
less on stability than would be estimated from the MO stability functions. 
Inclusion of the influence of subcanopy stability on subcanopy transfer is 
expected to improve formulation of transport of heat, momentum and pre- 
sumably trace gases such as carbon dioxide. 

Subcanopy drainage flows with strong capping inversions are common at 
the pine site and contribute to deviations from simple similarity theory. Ac- 
curate formulation of the within-canopy transport for these situations would 

23 


70 


require a separate subcanopy momentum equation. We are collecting data 
with better vertical resolution to determine the influence of this inversion on 
decoupling between surface fluxes and fluxes above the canopy. 

In spite of these difficulties and uncertainties with the above analyses, 
it is possible to improve upon existing simple formulations of within-canopy 
mixing, which completely neglect the dependence of within-canopy mixing on 
the stability. For example, an eddy diffusivity for heat has been constructed 
(Eq. 13) from the mixing length for neutral conditions (Eqs. 4, 5), the 
usual Monin-Obukhov stability functions and the formulation of the profile 
of the within-canopy friction velocity scale (Eqs. 14, 15), which depends on 
downward accumulated LAI and within-canopy stability. 

More complete subcanopy flux profiles and their microscale variability 
are needed to provide more dependable profiles of the mixing length for heat 
and the dependence of the nondimensional shear on canopy architecture and 
within-canopy stability. 
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Introduction 


This document describes the methods used to arrive at annual NEE, GEP 
and respiration estimates. General procedures common to all sites and all 
years are discussed in Sections 1-6. Site and year specific processing details 
follow the general discussion. 

This is a working document to describe the measurements and data pro- 
cessing that lead to the standard Ameriflux data set for the Metolius OR 
sites. It is not intended to be a manuscript for publication, nor does it 
include ongoing research efforts, although features of special interest along 
with some interpretation is presented in the discussion subsection of each 
site-year section. The biological measurements are not directly discussed in 
this document. 

The goal is construct quantitative and representative estimates for net 
ecosystem exchange of carbon (NEE) and to partition NEE into GEP (pho- 
tosynthesis) and respiration components. However, there are numerous po- 
tential problems. Long periods of missing data that coincide with a particular 
stage of biological activity could contribute a bias in the estimates. The flux 
instruments have quality problems in wet weather, which includes most of 
the winter season in the study area. The open-path LICOR is subject to 
calibration drift as the condition of the surface of the windows changes. 

The flux measurement height may be in the roughness sublayer, espe- 
cially in convective conditions, and there is no information on the horizontal 
or vertical variation of the flux. The single storage measurement site could be 
unrepresentative due to preferred locations of nocturnal carbon dioxide accu- 
mulation. There could be a significant mismatch between the flux footprint 
and the storage footprint, especially in stable conditions. The fetch changes 
with wind direction and stability, and the wind direction and stability vary 
systematically with season and time of day. Advection may be significant, 
but is not known. 
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Data collection 


Fast response wind and temperature data from a sonic anemometer (Camp- 
bell CSAT3) and fast response number density data from a colocated open- 
path C02/H20 analyser (LICOR-7500) positioned above the canopy are 
logged and saved at 10 or 20 hz for subsequent analysis. The mean C02 con- 
centration profile is measured at multiple levels using a single closed-path 
gas analyser (LICOR-6262) with multiple switching ports. Ancillary data 
can include mean temperature and wind profiles, radiation terms, soil prop- 
erties, precipitation and a micronetwork of meterological measurements in 
some cases. 


3 Quality control 

A complete description of this program is at 

http : //big . oce . orst . edu/Sof tware/QC_v3/guide/guide . html . 

The qc program is applied to raw time series of fast response variables used in 

eddy correlation flux calculations. This variable set includes ( u , v, w, T v , q, co2). 

The program performs the following steps: 

• Fast response specific humidity q ( gg -1 ) and co2 concentration (ppm) 
are computed from fast response number density (millimole m -3 ) mea- 
surements and fast response ambient air density (temperature and pres- 
sure fluctuations). This is sometimes called a point-by-point WPL 
“correction” . 

• Small segments of missing data in time series are replaced using linear 
interpolation. Small is defined as small enough such that the flux is 
not significantly effected, based on previous experience with other data 
sets. 

• The times series are despiked. The spike threshold is 3.25 times the 
local standard deviation. Continious spikes for periods longer than 5 
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seconds are considered real events and not spikes. Spikes are replaced 
using linear interpolation. Some data that is flagged by the qc test- 
ing is subjected to a “heavy-despiking” algorithm that removes longer 
sequences of spikes. 

• A series of tests are imposed that check the variable set for an un- 
physical range of values, unusally small or large variance, skewness and 
kurtosis, and for discontinuites in the mean. When suspect data is 
found it is masked out (set to missing data value). 

• A fast response “clean” data set is written for subsequent analysis. 

The LICOR-7500 has data quality problems when moisture (rain, snow, dew) 
or particulates accumulate on the windows. The CSAT generally has less 
problems than the LICOR-7500. 

Any quality control process is not 100% reliable. It is sometimes difficult 
to distinguish between instrument problems and rare but plausible physical 
behaviour. 


4 Second generation 

A complete description of this program is at 

http : //big. oce . orst . edu/Sof tware/2nd_and_3rdgen/ . 

This program computes means, variances, fluxes and spectra. Fluxes are 
computed for 3 different local averaging time scales; 100 seconds, 300 seconds 
and 600 seconds. For Ameriflux we use the 600-second fluxes. There is 
no detrending and fluxes are calculated using unweighted, non-overlapping, 
box-car averages. Fluctuations (e.g. w') for all points in the first 600-second 
window are computed by subtracting the 600-second mean. The same is 
done for co2. The product (e.g. w'co2') is then computed for all points in 
the window and then averaged over the window. This is repeated for the 
next 600-second window, etc. 
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No corrections are made for potential high frequency flux loss due to path- 
length averaging by the sonic or the LICOR, or for potential high frequency 
flux loss due to physical seperation between the sonic and the LICOR. 

Multiresolution cospectra and spectra are computed and saved for anal- 
ysis of the scale-dependence of the fluxes. 

We note that use of a 600-second local averaging time scale for computing 
fluxes in strongly stratified conditions may include poorly sampled mesoscale 
motions. As a result, the flux is characterized by large random sampling 
error. The contribution to the calculated flux at these large time scales (low 
frequencies) may be larger than the contribution from turbulent scales in 
very weak turbulence. 

4.1 Tilt correction 

A tilt correction is applied to the wind components prior to computing fluxes. 
The tilt correction can be time-dependent to accomodate changes in the sonic 
orientation due to any physical realignments that may take place during the 
year. In this case, the entire procedure (see below) is repeated seperately for 
each period. Individual hours with weak horizontal winds (less than 2 ms -1 ) 
are excluded from the tilt correction development procedures since the tilt 
angle can be erratic in weak winds. 

A practice constant offset is removed from the vertical motion making w 
zero for the entire period. The removal of the offset does not directly affect 
the computed flux, since only a constant is removed for each record. However, 
we note that applying a subsequent tilt rotation to the data with no prior 
removal of the offset, converts the vertical motion offset into horizontal flow. 

After removing the practice offset, a practice tilt angle is computed for 
each 1-hour record which eliminates the record mean vertical motion. These 
tilt angles are averaged over the entire period for a sequence of wind direction 
categories, giving wind direction-dependent average tilt angles. When the 
practice tilt angles appear to be consistant with a tilted anemometer, the 
actual tilt correction is applied. For each 1-hour data record, a lookup table 
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is used to determine the offset and tilt angle to use depending on the time 
and the wind direction, and the rotation is done on the fast response wind 
components. The process includes a horizontal rotation to along-wind ( u ) 
and cross-wind ( v ) components such that u > 0 and v = 0 to distinguish 
between along- and cross-wind components of the wind stress. This method 
does not force the record-mean vertical motion to zero. The correction is 
applied to all records regardless of wind speed. 

5 Third generation 

A complete description of this program is at 

http : //big . oce . orst . edu/Sof tware/2nd_and_3rdgen/ . 

This program calculates quantities related to evaluating Monin Obukov sim- 
ilarity theory. Several estimates for different types of flux sampling errors 
are computed including random sampling error, systematic error, a flux in- 
termittency parameter and flux non-stationarity. 

Half-hour mean quantities axe computed from the 2nd generation 600- 
second data. Three 10-minute average values are averaged to make one 30- 
minute average. The sensible heat flux is calculated from the sonic heat flux 
(which is close to the virtual temperature flux) and the moisture flux. 

6 Ameriflux processing 

6.1 Flow distortion 

The fluxes for half-hourly periods where the flow is potentially disturbed due 
to sonic support structures, other instruments or the tower itself are masked 
out (set to missing data value). The width of the wind sector deemed to have 
potential flow distortion is specific to each site and depends on the type of 
tower, boom length and configuration. 
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6.2 Outliers 


Flux outliers are removed to avoid their influence on gap-filling. Values of 
the eddy correlation (EC) flux outside a user specified range are masked out. 
The range is selected after consulting a time series plot of all the fluxes, 
and should include only a very small fraction of the data. It is likely that 
some of these outliers are associated with instrument problems that were not 
detected by the quality control testing. 

6.3 Gap-filling 

The steps in the gap-filling procedure are outlined below. 

• 1. Temporal fill wide gaps using linear interpolation. A wide gap is 
a sequence of consecutive days with a very high percentage of missing 
data (90% or more). Typically, the percent of missing data in these 
cases is 100%. When equivalent width data segments just prior to and 
subsequent to the wide gap segment have a lower percent missing data 
(60% or less), then linear interpolation is used to fill across the wide gap. 
Means are computed for the prior and subsequent segments for each 
of 48 half-hour periods. The half-hourly values are then interpolated 
across the wide gap. This procedure is implemented using a sequence 
of block windows of width 35,30,-10,5 days. 

• 2. Vertical fill 1-pt and 2-pt gaps using linear interpolation for co2(z) 
and T(z). A 1-pt gap is one vertical level with missing data with good 
data both above and below, for the same time period. 

• 3. Temporal fill 1-pt and 2-pt gaps using linear interpolation for all 
variables. A 1-pt gap here is one half-hour with missing data with 
good data for the prior and subsequent half-hour periods. 

• 4. The mean diurnal averages are calculated for all variables using 
all non-missing data within a 15 day window width. The mean diurnal 
average consists of the averages over the window for each of 48 half-hour 
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periods during the day. Missing data within the window is replaced 
with the correpsonding average over the window if at least 50% of the 
data for that period were included in the average. For example, if a 
particular half-hour period is missing for 8 or more of the 15 days in the 
window, the missing data for that half-hour period will not be replaced 
by the average. 

This entire process is then repeated partitioning the data into 8 three- 
hour periods instead of 48 half-hour periods. The probability of finding 
enough data to meet the 50% criteria is greatly increased because the 
time-of-day criteria is relaxed from a half-hour to three-hours. In gen- 
eral, this will depend on the character of the missing data and could 
be site-specific. When replacing missing data with the average, the six 
half-hour averages in the three-hour period are all assigned the same 
value. 

• 5. Repeat step 4 for increasing window sizes of 30, 45, etc. days until 
all data is gap-filled for all variables. 

6.4 Storage term 

The storage of co2 and heat in the air between the ground and the flux 
measurement height is computed by vertically integrating the local time ten- 
dancy. The time change is computed using centered (2 St) differences and 
half-hour means. The vertical integration is done assuming a piecewise linear 
fit to the co2(z) and T(z) profile data. The profiles are extrapolated from 
the first measurement level down to the surface using the slope between the 
first and second levels. 

6.5 w*-filtering 

Nocturnal measurements of the co2 EC flux plus storage sometimes increase 
with increasing friction velocity (u„). Since respiration is not thought to 
be function of «*, this indicates that other processes may be acting. The 
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generally accepted interpretation is that horizontal or vertical advection of 
co2 must be important. Another possible interpretation is that the storage 
measurement is unrepresentative, possibly due to horizontal heterogentiy of 
the storage. Another is that the flux footprint region could be quite different 
than the footprint of the storage measurement in stable flows. 

One approach to this problem is to apply filtering, where a model 
of respiration is developed using the nocturnal F c + S c measurements in 
only strong mixing (high «*) conditions. In such conditions, advection and 
horizontal heterogentiy are probably less important. The critical u, value is 
site dependent. The friction velocity is calculated as the square root of the 
magnitude of the 30-minute averaged wind stress components 

u t = ( w'u ' 2 + w'v ,2 ) 1 / 4 . (1) 

The critical u t value is estimated in off-line mode by exmaining the in- 
dependence of F c +S c . The year is broken up into a sequence of 60-day periods 
to accomodate seasonal differences. Plots of nocturnal F c + S c bin-averaged 
by u* category are examined to see if and where a levelling off occurs with 
increasing it*. Such plots are included for each site in the latter part of this 
document. Ideally, F c + S c should not be sensitive to it* for it, larger than 
the critical value. We define “nocturnal” as the half-hour periods where the 
solar zenith angle exceeds ninety degrees. 

Sensitivity tests indicate that the F c + S c vs (it*) relationship is heavily 
influenced by outliers. A more robust relationship is found after discarding 
points in the upper 2% and lower 2% of the frequency distribution of both 
F c + S c and it*. Negative nocturnal values of F c + S c are discarded. Such 
values are suspect probably due to large random flux sampling errors which 
are more severe in stable conditions. 

Once the critical it, is determined, nocturnal F c +S c for high it, conditions 
is related to predictors of respiration, such as subcanopy air temperature, soil 
temperature and soil moisture, to develop a model of respiration. The precise 
form of the model depends on the special characteristics of the site and the 
availablity of measurements. When available, soil chamber measurements of 


9 


84 


respiration are consulted in the model development. Specific details of the 
respiration model are discussed below in the individual site-year sections. 

A potential problem with the “tt ,,-filter respiration modelling” approach 
is that the warmest soil temperatures occur in summer during the day with 
weak mixing (low u„) conditions. The respiration during such conditions 
can not be estimated using F c + S c , and thus these conditions can not be 
included in developing the model coefficients. To estimate respiration for such 
conditions, the model must be used to extrapolate to higher temperatures, 
although the reliability of such an extrapolation is not known. 

Respiration is estimated by applying the model to all daytime and noctur- 
nal periods regardless of mixing strength. Recall that the model coefficients 
were derived using the F c + S c measurements for only the high u* nocturnal 
periods, and can be specified seperately for each 60-day period. 

GEP is calculated as F c + S c minus the modelled respiration during the 
day, and is set to zero at night. NEE is then given by GEP + respiration. 
Respiration, GEP and NEE are defined as negative for carbon uptake by the 
ecosystem and positive for release of carbon to the atmosphere. 

Annual NEE ( gCm ~ 2 ) is computed as 

NEE = N d 86400 10“ 3 (12/44) A (2) 

where A is the annual average NEE (mg co2 m~ 2 s~ x ) and N d is the number 
of days per year. The same approach is used for computing annual sums of 
respiration and GEP. 

7 IP03 

The first application of the OSU-BLG processing was for the 2003 data from 
the intermediate pine site near Sisters, OR, USA. The above canopy flux mea- 
surements were made 31 m above ground using a CSAT3 sonic anemometer 
and open-path LICOR-7500. The average height of the canopy was 17 m, so 
the flux instruments were approximately 14 m above canopy, or at 1.8 canopy 
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heights. Significant solar radiation reaches the forest floor underneath the 
canopy in summer indicating a fairly sparse canopy. During summer nights, 
the subcanopy layer remains significantly stratified even with the highest 
above canopy u* conditions. 

Mean co2(z) (closed path LICOR-6262) was measured at 1, 4 and 30 m. 
Mean temperature (HOBOs) was measured at 3, 6, 10, 20 and 30 m. Mean 
horizontal winds (Handars) were measured at 3, 6, 10, 20 and 30 m. The 
soil heat flux was measured at 5 locations and averaged. Soil temperature 
was measured at 6 depths between 2 and 64 cm at one location. Soil mois- 
ture content was measured at one location. Automated daily mean chamber 
measurements of soil respiration were made during April through Spetember. 

7.1 Data coverage 

A total of 6790 hours of data were logged (78% annual coverage). Approx- 
imately 10% of these hours contained enough missing data that they could 
not be saved and were discarded. Approximately 10% of the surviving hours 
were discarded by quality control tests. Problems with either the co2 or h2o 
measurement (usually both) from the LICOR-7500 were 25 times more likely 
to occur than problems with the CSAT. 

Subsequent to quality control, a problem with the 31-m fluxes was dis- 
covered for the period 11-Sep to 29-Sep (DOY 255-272), and these data were 
discarded. An additional problem with the 31-m fluxes was discovered dur- 
ing the tilt correction phase of the processing for the period from 4-Nov to 
the end of the year, and these data were discarded. Apparently, the mount 
holding the sonic anemometer to the boom came loose. 

The gap-filling did not work very well for the sensible and latent heat 
fluxes and the co2 flux during November and December, where there was 
almost no good quality flux data. The large amount of missing data (Table 
1) sever ly tested the gap-filling procedures. 


11 


Table 1. IP03 percent of available EC co2 flux (F c ) and storage (S c ) data. 


Month 

J 

F 

M 

A 

M 

J 

J 

A 

S 

O 

N 

D 

F c 

38 

68 

40 

46 

21 

74 

91 

81 

31 

76 

5 

0 

S c 

88 

100 

98 

99 

90 

33 

72 

0 

83 

56 

54 

42 


7.2 Respiration modelling 

Nocturnal F c + S c increases with increasing u* at this site during some pe- 
riods of the year (Figure 1). Based on these data, we select a critical u » 
value of 0.5 ms~ l for the purpose of applying the ^.-filtering, and apply this 
constant value for the entire year. The respiration from nocturnal EC flux 
plus storeage measurements during high u * periods agrees reasonably well 
with the automated chamber measurements (horizontal lines in Figure 1). 
The eddy correlation estimates of respiration are slightly higher than the 
chamber values presumably due to respiration from the foilage, which would 
not be captured by the chamber measurements. 

The subcanopy air temperature, 2-cm soil temperature and 4-cm soil tem- 
perature were all reasonably good predictors of respiration (measured noc- 
turnal F c + S c in high u* conditions). We choose the 4-cm soil temperature. 
The temperature-dependence was fit to the form 

Res = ot T? (3) 

using least squares regression. The respiration was first bin-averaged by 
temperature category (circles in Figure 2) and then the regression was done 
on the bin-averaged data to obtain the coefficients, a = 0.037 and 0 = 0.563 
and the model (solid line Figure 2). The model coefficients change slightly to 
0.033 and 0.595, respectively, using a critical u* of 0.7 ms -1 instead of 0.5. 

Theoretically, the respiration from the soil should increase exponentially 
with soil temperature and decrease with decreasing soil moisture content. 
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The reason that soil temperature alone describes the observed respiration so 
well at the Metolius sites is that soil moisture content is negatively correlated 
with soil temperature, at least for the first half of the year. During July the 
soil moisture content reaches a minimum, and then does not decrease any 
further through the late summer and fall. The implicit inclusion of soil 
moisture content in the soil temperature leads to the relatively weak soil 
temperature dependence of the respiration found here. 

The daily average, temperature-dependent respiration model agrees quite 
well with the daily average, corrected, automated soil chamber measurements 
through the end of June, but then is significantly larger than the chamber 
measurements during July through September. Adding soil moisture content 
to the respiration model does not significantly change this result. This could 
be due to the soil moisture measurement “pegging” at a low value during 
July and remaining relatively unchanged until the late fall rains begin. The 
soil moisture is measured using the CS615, which when installed vertically, 
gives a measure of the soil moisture content in the upper 30 cm of soil. 
Another possibility is that the fractional respiration coming from the foilage 
is higher in late summer than in early summer. The chamber measurements 
miss foilage respiration. 

As a sensitivity test, the respiration model coefficients were calculated 
seperately for each of six 60-day window time periods during the year. This 
approach allows a different temperature dependence for different seasons of 
the year. The motivation is that while soil temperatures are similar in early 
and late summer, the soil moisture content is significantly less in late sum- 
mer. The comparison with the automated chamber data is very similar to 
that described above for the simpler respiration model with constant coeffi- 
cients for the whole year. That is, the model respiration is in good agreement 
with the automated chamber measurements through June, but then is larger 
than the chamber measurements for July through September. The respira- 
tion during late summer and fall is slightly reduced for the seasonal model 
compared to the constant coefficient model, and is therefore slightly closer 
to the chamber measurements. Since the approach using constant year-long 
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coefficients appears to be just as reasonable as the more complex model with 
seasonal-dependent coefficients, we choose to apply the simpler one. 

Research is ongoing on the respiration. 

7.3 Annual estimates 

The estimate of annual NEE at this site is perhaps high compared to ex- 
pectations for a semi-axid ponderosa pine site (Table 2). As a sensitivity 
test, increasing the u* threshold from 0.5 to 0.7 ms -1 slightly decreases the 
respiration, and thus increases the annual carbon uptake, contrary to ex- 
pectations. The fact that increasing the critical u * from 0.5 to 0.7 does not 
increase the respiration is not conclusive since: 1) the amount of high u r data 
is limited and 2) the footprint increases with wind speed. 

As another sensitivity test, discarding the u*-filter approach and using 
the automated chamber measurements for respiration, decreases the annual 
respiration and increases the annual carbon uptake. In this calculation, the 
respiration model with critical u » = 0.5 ms -1 is used when the chamber data 
is missing. The chamber measurements are not extrapolated into the winter 
season. 

The July mean diurnal cycle of the terms in Table 2 is shown in Figure 
3, while the annual cycle of weekly-averages is in Figure 4. These fields are 
based on the standard calculation with w*-filtering and a threshold value of 
0.5 ms -1 . 


Table 2. IP03 annual estimates ( gCm 2 ). 


F c 

So 

GEP 

Re. 

NEE 

-513 

0 

-1302 

936 

-366 
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7.4 Discussion 


An under-estimate of respiration would contribute to the large negative NEE 
indicating relatively large carbon uptake compared to expectations for this 
site. Comparisons with the automated chamber measurements are incon- 
clusive. In June and July with the warmest soil temperatures, the cham- 
ber measurements of respiration exceed our modelled estimate, however, the 
chamber estimates are significantly lower than the model estimates in the 
early and late summer. It is not clear why the chamber measurements drop 
off so sharply in late July and early August, while the soil temperature and 
soil moisture are not changing rapidly. Spatial heterogenity is an issue with 
the chamber measurements. 

A possible reason for the relatively large NEE estimate is that even with 
strong mixing above the canopy (high u t conditions), the mixing underneath 
the canopy is still suppressed due to a suprisingly strong subcanopy temper- 
ature stratification. The subcanopy stratification decreases with increasing 
w» above the canopy as expected, but does not go to zero. In summer, the 
stratified layer results from strong radiative cooling at the surface due to 
clear skies, dry soil and a sparse canopy. Because of the suppressed mixing 
in the subcanopy layer, the nocturnal flux measurements made above the 
canopy may under-estimate the respiration even under the strongest mixing 
conditions observed. 

Another potential complication is that when the wind is from the south 
or southwest, it comes from a region with smaller LAI. The tower site is 
in an area of local maximum LAI. At night, the surface footprint is larger 
(enormous with subcanopy stratification) and is expected to be smaller with 
daytime heating. Therefore, the daytime carbon uptake is dominated by the 
high LAI area in the vicinity of the tower, while the nighttime respiration is 
more strongly influenced by areas with lower LAI. 
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8 IP04 


The second application of the OSU-BLG processing was for the 2004 data 
from the intermediate pine site near Sisters, OR, USA. See the IP03 section 
for details. Differences between the two years are discussed here. 

More levels of mean co2 and temperature were deployed on the flux tower 
in 2004. Mean co2(z) (closed path LICOR-6262) was measured at 1, 3, 6, 15 
and 31 m, and mean temperature (HOBOs) was measured at 3, 6, 10, 15, 20 
and 30 m. 

8.1 Data coverage 

2004 has less missing data compared to 2003 (Table 3). For the entire year, 
a valid co2 EC flux measurement was logged and passed the quality control 
testing for 82% of the 30-minute periods. The correpsonding number for the 
co2 storeage measurements is 88%, with nearly all the missing storeage data 
being in January. There is no flux data prior to 18-January. 


Table 3. IP04 percent of available EC co2 flux (F c ) and storage ( S c ) data. 


Month 

J 

F 

M 

A 

M 

J 

J 

A 

S 

O 

N 

D 

F c 

25 

74 

79 

76 

84 

83 

90 

83 

86 

83 

76 

67 

Sc 

3 

92 

99 

95 

98 

72 

99 

100 

100 

99 

99 

100 


8.2 Respiration modelling 

Based on the ^.-dependence of F c + S c , we choose a critical u, threshold of 
0.55 ms -1 and apply this constant value for the entire year (Figure 5). The 
^.-dependence is similar to what was found at this site in the previous year. 

The temperature-dependence of the respiration (measured nocturnal F c + 
S c in high it* conditions) was fit to 
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R es = a T p 


( 4 ) 


using least squares regression where T is the 4-cm soil temperature. 

The model coefficients (a = 0.054 and 0 — 0.489) tend to over-estimate 
respiration for cooler soil temperatures, and under-estimate respiration for 
warmer temperatures (Figure 6). The model fit is sensitive to the warmest 
soil temperatures included in the regression, here about 25 C. The data 
suggest that the respiration actually decreases or levels off with increasing soil 
temperature for soil temperatures above 20 C. This decrease is unexpected 
based on our understanding of soil respiration, and could be due to random 
sampling problems caused by an insufficient data sample size. The general 
problem discussed in Section 7.2 also applys here. 

8.3 Annual estimates 


Table 4. IP04 annual estimates (gCm 2 ). 


F c 

5c 

GEP 

Res 

NEE 

-604 

0 

-1479 

1142 

-337 


The July mean diurnal cycle is shown in Figure 7, while the annual cycle 
of weekly-averages is in Figure 8. 

8.4 Discussion 

The annual NEE for 2004 (337 gCrrT 2 ) is about 8% less than that found for 
2003 (366 gCm ~ 2 ). Both GEP and respiration increased in 2004 compared 
to the previous year. 
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9 YP04 


This is the new young pine site near Sisters, OR, USA. The above canopy flux 
measurements were made 12 m above ground using a CSAT3 sonic anemome- 
ter and open-path LICOR-7500. The average height of the canopy was 3 m, 
so the flux instruments were at 4 canopy heights. The short pines are uni- 
formly and sparsely distributed and the canopy is not closed. 

Mean co2(z) (closed path LICOR-6262) was measured at 1, 5, 12 and 
18 m. The 18 m mean co2 measurements did not begin until June. Mean 
temperature (HOBOs) was measured at 1, 2.5 and 5 m. The soil heat flux 
data is bad until the very end of the year, when two new heat flux plates 
were installed. Soil temperature was measured at 3 depths between 2 and 
32 cm at one location. Soil moisture content was measured at one location. 
There is no pyranometer at this site. 

9.1 Data coverage 

The missing EC flux data is fairly evenly distributed through the year, and 
there are no months with more that 50% missing data (Table 5). This type 
of missing data is better suited to the gap-filling methods. There is very 
little missing co2 profile data. 


Table 5. YP04 percent of available EC co2 flux (F c ) and storage (5 C ) data. 


Month 

J 

F 

M 

A 

M 

J 

J 

A 

s 

O 

N 

D 

F c 

51 

66 

77 

80 

80 

89 

66 

77 

78 

74 

77 

59 

S c 

97 

100 

100 

98 

91 

89 

99 

100 

99 

99 

99 

100 
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9.2 Respiration modelling 

There is no clear ^.-dependence of F c + S c at this site, presumably due to 
the low LAI and sparsely distributed short pines (Figure 9). The shortness 
of the trees also helps eddies communicate with the surface. That is, there 
appears to be strong coupling between the 12-m measurement level and the 
surface, even on weak mixing nights. 

The temperature-dependence of the respiration (measured nocturnal F c + 
S c ) was fit to 

Re, — a T& (5) 

using least squares regression where T is the 2-cm soil temperature. (Figure 
10). The coefficients are a — 0.044 and (3 = 0.410 based on the 2-cm soil 
temperature (there are no 4-cm soil temperature measurements). Because 
there is no critical u* value, all the nocturnal data were used to develop the 
respiration model coefficients, not just the high u t nocturnal data. 

9.3 Annual estimates 


Table 4. YP04 annual estimates ( gCm 2 ). 


F c 

5c 

GEP 

Re, 

NEE 

-233 

0 

-981 

878 

-103 


Annual NEE at the young site is about one-third of that found for the 
intermediate site. The July mean diurnal cycle is shown in Figure 11, while 
the annual cycle of weekly-averages is in Figure 12. 
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9.4 Discussion 

In contrast to the intermediate site, there is no indication that respiration is 
“missed” by the EC flux measurement at 12 m for low u, nocturnal conditions 
at the young site. 

There is evidence of early morning vertical flushing of co2 associated with 
upward EC co2 flux shortly after sunrise. This is accompanied by a sharp 
decrease in mean co2 concentration. This also coincides with a shift in wind 
direction from southwesterly to northeasterly, which could be associated with 
a thermally induced upslope flow. The terrain slopes down to the east and 
northeast. 


10 Intercomparisons 

The annual course of the cummulative NEE ( gCm ~ 2 ) for each site-year is 
shown in Figure 13. The cummulative NEE is computed as 

Cummulative NEE = 1800 10~ 3 (12/44) NEEac (6) 

where NEE ac is the cummulative sum of variable NEE in the data file, where 
NEE in the data file has units of (mg co2 m~ 2 s~ l ) and represents a 30-minute 
time average. 
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Figure 1: IP03 seasonal pattern of bin-averaged nocturnal F c + S c ( mg co2 
m -2 s -1 ) versus friction velocity (ms -1 ). Solid horizontal lines are estimates 
of daily mean soil co2 efflux (respiration) from automated chamber mea- 
surements, corrected using a relationship between manual and automated 
chamber data. 
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Figure 5: IP04 seasonal pattern of bin-averaged nocturnal F c + S c (mg co2 
m -2 s -1 ) versus friction velocity (ms -1 ). Solid horizontal lines are estimates 
of daily mean soil co2 efflux (respiration) from automated chamber mea- 
surements, corrected using a relationship between manual and automated 
chamber data. 
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Figure 9: YP04 seasonal pattern of bin-averaged nocturnal F c + S c (mg co2 
m _2 s _1 ) versus friction velocity (ms '). 
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Abstract 


Numerous recent studies calculate horizontal and vertical advection terms 
for budget studies of net ecosystem exchange of carbon. One potential un- 
certainty in such studies is the estimate of mean vertical motion. This work 
addresses the reliability of vertical advection estimates by contrasting the 
vertical motion obtained from the standard practise of measuring the verti- 
cal velocity and applying a tilt correction, to the vertical motion calculated 
from measurements of the horizontal divergence of the flow using a network 
of towers. Results are compared for three different tilt correction methods. 

Estimates of mean vertical motion are sensitive to the choice of tilt cor- 
rection method. The short-term mean (10 to 60 minutes) vertical motion 
based on the horizontal divergence is more realistic compared to the esti- 
mates derived from the standard practise. The divergence shows long-term 
mean (days to months) sinking motion at the site, apparently due to the 
surface roughness change. Because all the tilt correction methods rely on the 
assumption that the long-term mean vertical motion is zero for a given wind 
direction, they fail to reproduce the vertical motion based on the divergence. 


Keywords 

Vertical Advection, Mean Vertical Motion, Tilt Correction, Horizontal 
Divergence 
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1 Introduction 


Early studies estimating net ecosystem exchange of carbon (NEE) typically 
included only the storeage and eddy flux components. It is now generally 
accepted that this method leads to what appears to be an underestimate of 
nocturnal respiration of carbon dioxide on nights with weak mixing. One 
possible interpretation of the missing carbon dioxide is that horizontal or 
vertical advection is important (Lee, 1998; Baldocchi et al., 2000; Aubinet 
et al., 2003; Staebler and Fitzjarrald, 2004). An empirical patch for this 
problem is the so called u, filter approach, where the measured eddy flux 
during weak turbulence nocturnal periods, where advection is most likely to 
be important, is replaced with a model of the respiration based on the eddy 
flux during strong mixing nocturnal periods (Pattey et al., 2002). 

The method of calculating NEE proposed by Lee (1998) added a vertical 
advection term to the storeage and eddy flux terms in the NEE budget. This 
approach was applied by Baldocchi et al. (2000) who reported that including 
vertical advection improved their budget closure, however, they inferred that 
horizontal advection may also be important. Aubinet et al. (2003) concluded 
that horizontal advection of carbon dioxide was mostly cancelled by verti- 
cal advection. Their hypothesis was that horizontal and vertical advection 
are linked by an entrainment mechanism where the air above the canopy is 
brought downward into the subcanopy in drainage flows. This suggests that 
it is inappropriate to include only the vertical component of advection in 
the budget. Staebler and Fitzjarrald (2004) found that vertical advection of 
carbon dioxide was small while horizontal advection was significant. Includ- 
ing horizontal advection improved their NEE budget closure, especially in 
the summer. Feigenwinter et al. (2004) found that horizontal and vertical 
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advection of carbon dioxide tended to cancel at night, yet pointed out that 
the large scatter in the advective fluxes needs further investigation. 

All of the above studies estimate mean vertical motion by applying a tilt 
correction (Section 3) to the vertical velocity measured by a 3-dimensional 
sonic anemometer, however, it has not been demonstrated that such an ap- 
proach is capable of resolving a small time- averaged vertical motion. Can the 
true vertical motion be extracted from the measurements which are contami- 
nated by flow distortion, tilt of the sensor and sloping terrain? In this study, 
we investiagate the feasibility of quantifying mean vertical motion, and thus 
vertical advection, using sonic anemometer measurements. 

A network of towers were deployed to provide two independent estimates 
of mean vertical motion; 1) from measurement of vertical velocity at a cen- 
tral tower, and 2) from the horizontal divergence using a network of towers 
surrounding the central tower. Prom incompressible mass continuity, the 
time-averaged vertical velocity based on the divergence is given by 


f z=h ,du dv. , 

w{h)= -L { di + ^ )dz (1) 

using the usual notation where overbars denote a time-average and «J(0) =0. 
Multiple measurement levels of the horizontal flow on the surrounding tow- 
ers provide some details on the vertical structure of the horizontal gradients 
of the mean wind. Two levels of vertical velocity measurements above the 
canopy on the central tower provide a consistancy check and allow identifica- 
tion of complex flow situations where the magnitude of the vertical velocity 
may not increase with height. 

A strong relationship between the mean vertical motion from the two 
independent methods would provide confidence for estimates of vertical mo- 
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tion, and subsequent estimates of vertical advection. However, when the 
two estimates substantially disagree, such confidence can not be established. 
This approach is not definitive due to potential errors in both estimates of 
w. The comparison is complicated by the fact that w at a point responds to 
divgerence on a large range of spatial scales while our network for calculating 
the divergence resolves only one spatial scale. We attempt to address this 
problem by examining the relationships between the w estimates for a wide 
range of averaging times. 

2 Field experiment 

The measurements are from a semi-arid young ponderosa pine site in central 
OR, USA, during August, September and October of 2004. The young 3- 
m tall pine trees were uniformly and sparsely distributed across the site. 
Beneath the pines was bare sandy soil with sparse clumps of grass. A stand 
of much taller pines surround the site. The site was chosen in part due 
to the large roughness change and the expected acceleration of the mean 
flow (divergence) over the shorter trees. The boundary between the young 
and taller pines is roughly denoted by the edge of the area in Figure 1. 
Within the boundaries of the site, the vegetation could be classified as very 
homogeneous. The terrain slopes steadily upward to the southwest and west 
with 100 m elevation gain over the first 5 km, or a slope of 2%, and slopes 
weakly downward (< 1%) or is flat in other directions. More than one-half 
the time the wind direction was in the sector 230 to 300 deg (westerly winds). 

The central tower was equipped with 3-dimensional CSAT3 sonic anemome- 
ters (Campbell Scientific) at 5 and 12 m above ground recording the three 
wind components and the sonic temperature at 10 hz. The orientation of 


4 


114 


the anemometers was carefully measured with a SmartTool (MD Buliding 
Products), a digital instrument with a precision of 0.1 deg. The pitch, roll 
and azimuth of the 12-m anemometer were measured as 0.9, 0.0 and 272 
deg, respectively. Two-dimensional sonic anemometers (Vaisala Inc, formerly 
Handar) were deployed at 1, 5 and 12 m above ground on towers A and C 
and at 1 and 5 m on towers D, E and B for measuring the mean horizontal 
wind. All periods with flow from 60 through 150 deg were discarded due to 
possible flow distortion due to flow through the tower and support structures 
prior to reaching the instruments. 

The site is characterized by sustained synoptic high pressure and weak 
mean flow in summer and fall. For the 3 month experiment, the average wind 
speed was 0.7, 1.4 and 1.9 ms -1 at the 1, 5 and 12-m levels on tower B, re- 
spectively. The weak mean flow, clear skies and sparse canopy lead to strong 
radiational cooling of the ground at night, resulting in strong stratitifcation 
and very weak turbulence near the surface. 

Direct comparison between the mean wind speed measured by a CSAT3 
and a Handar is afforded by the dual deployment at 5 m on the central 
tower. The 10-minute vector average wind speed from the CSAT3 (Handar) 
was 1.36 (1.34) ms -1 for the entire experiment, and the two wind speeds 
were correlated at r = 0.98. 


3 Tilt correction 

3.1 Background 

A tilt correction is applied to 3-dimensional sonic anemometer measurements 
to remove; 1) apparent vertical motion due to tilt of the sensor from true 
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(gravitational) vertical, and 2) apparent vertical motion due to sloping topog- 
raphy. Both of these influences lead to an artifical wind direction dependence 
of the vertical velocity as measured by the anemometer. As noted by Paw 
U et al. (2000) and others, it is not possible to distinguish between these 
two effects from the vertical velocity measurements alone. A tilted sensor 
over flat terrain measures an apparent upward vertical motion for horizontal 
flow in one direction and apparent downward motion for flow in the opposite 
direction. The apparent vertical motion signature is identical for a perfectly 
levelled sensor over idealized sloping terrain with terrain following flow. Di- 
rectionally dependent flow distortion effects are typically assumed to be ac- 
counted for by the internal software in the anemometer, although Hogstrom 
and Smedman (2004) recently pointed out that the manufacturer’s correc- 
tions, which are based on low turbulence level wind tunnel studies, may not 
be optimum for the higher turbulence conditions found in the field. 

In earlier days of micrometeorological measurements, standard procedure 
was to rotate the three wind components into streamline coordinates for 
each 30-minute or 1-h time period individually (e.g., Kaimal and Finnigan, 
1994). This approach defines the short-term mean vertical motion to be zero. 
The main motivation for this approach was that the mean vertical motion 
measurements could not be trusted due to contamination by sensor tilt and 
elevation slope, limited sample size and possibly flow distortion. 

More recently, Lee (1998) and others identified the importance of non- 
zero w for calculating vertical advection of carbon dioxide. Even very small 
values of w can lead to significant advection in the presence of large vertical 
gradients, such as typically found near the surface for carbon dioxide on 
nights with weak mixing. Since the streamline coordinates approach exactly 
removes all short-term mean vertical motion, other methods were adopted 
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which allow non-zero short-term w. 

In the new standard approach, the long-term w m (<j>) averaged over weeks 
to months is removed from the short-term averages of w m , where <f> is wind 
direction and w m is the measured time-averaged vertical velocity (Lee, 1998; 
Baldocchi et ah, 2000; Paw U et ah, 2000; Wilzcak et ah, 2000; Finnigan 
et ah, 2003). The procedure is similar in concept to applying a high-pass 
filter to the mean vertical motion seperately for each wind direction. The 
approach assumes that the long-term mean vertical motion for a given wind 
direction is zero. The short-term deviations in w that remain after removing 
the long-term are presumed to correspond to short-term events of 

either horizontal convergence or divergence. Long-term measurements are 
required to implement this approach in order to obtain robust estimates of 
the long-term 

While the philosophy of removing the directionally dependent, long-term 
mean vertical motion is present in all current tilt correction methods, the 
details of applying it vary between studies. Lee (1998) proposed linear re- 
gression of w m on the mean horizontal wind speed seperately for each wind 
direction category, which accounts for the fact that the apparent vertical mo- 
tion due to sensor tilt is proportional to the horizontal wind speed. In the 
planar fit technique (Wilczak et al. ,2001), multiple linear regression of w m 
on the two horizontal wind components results in a tilted plane. This method 
applies to the idealized case where the surface is uniformly tilted (planar) 
with respect to the sensor, or equivalently, where the sensor is tilted over per- 
fectly flat terrain. A third approach examines the wind direction dependence 
of the tilt angle (Paw U et al., 2000; Vickers and Mahrt, 2003; Feigenwinter 
et al., 2004). The ^-dependence of the long-term average tilt angle can be 
represented either by the bin-averaged estimate or using a regression fit to 


7 


117 



azimuthal angle. When the (/(-dependence can be approximated by a simple 
sinusoidal function of azimuth angle, this method is identical conceptually 
to the planar fit technique (Paw U et al., 2000). In this sense, the planar fit 
approach is a special case of the tilt angle method. For non-idealized terrain, 
a planar fit may not be appropriate, and the bin-average (/(-dependence can 
be used. 

The regression used by Lee’s approach and the planar fit method exactly 
remove the mean vertical motion averaged over the entire experiment, while 
the tilt angle method does not. Lee’s approach removes the long-term w for 
each wind direction category individually, while the planar fit and tilt angle 
methods do not. 

3.2 Problems 

As pointed out by Lee (1998), a potential problem is that the tilt correction 
may remove real vertical motion associated with local circulations driven by 
surface heterogenity. Examples include a change in surface roughness, dif- 
ferential daytime heating and systematic diverging nocturnal drainage flows. 
The site here is influenced by a surface roughness change, however, it is not 
known if systematic circulations induced by daytime differential heating are 
important. It is also not known if the drainage flow is divergent (accelerat- 
ing) or convergent (decelerating) at the site. For westerly flow, the site is 
approximately at the bottom of a slope in that the elevation increases faster 
to the west than it decreases to the east. A complication with identifying 
the drainage flow at this site is that the predominant flow is from the same 
direction as the drainage flow. 

Given the rough-to-smooth surface roughness change for all wind direc- 
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tions, and the persistant synoptic high pressure pattern, we would expect 
overall mean sinking motion at the site. Such long-term mean vertical mo- 
tion would be at least partially removed by any of the tilt correction methods 
described above. 

We note that drainage flow is not a sufficient condition for horizontal 
divergence and mean sinking motion. Drainage flows need to converge some- 
where. This is an important point because nearly every NEE study in the 
literature reports mean sinking motion at night and attributes it to drainage 
flows. If this were indeed the case, it indicates that the site selection criteria 
were biased in that all sites are in a similar nocturnal flow regime. Mean 
sinking motion coupled with decreasing carbon dioxide concentration with 
height results in a vertical advection term that is the correct sign to explain 
the missing carbon dioxide found at most sites. However, because other 
terms in the NEE budget, such as horizontal advection and horizontal flux 
divergence, are typically not known, it is impossible to make definitive con- 
clusions. The typical reported finding of rising motion during the day and 
sinking motion at night could be an artifact of the tilt correction method. 

3.3 Methodologies 

We apply three different tilt correction methods to the data. The same pro- 
cedures are applied independently to the 5 and 12-m measurements although 
we focus on the 12-m estimates of w. All time mean quantities in this sec- 
tion, denoted with an overbar, represent 100-s time averages. The choice of 
100-s averaging includes wind direction variability on scales down to 100 s. 
Longer averaging times would miss additional meandering of the mean wind 
which could contaminate the wind direction dependence of w. It is difficult 
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to say, because many studies do not include the details, but it appears that 
studies in the NEE literature often employ a 30-minute averaging time to 
define the mean wind direction for the purpose of developing the tilt cor- 
rection. Berger et al. (2001) use the 1-h average wind direction. Wilczak 
et al. (2001) propose using 5-minute average wind components to derive 
the regression coefficients in the planar fit approach. We note that the tilt 
correction parameters could be developed using the high rate (e.g. 10 hz) 
data, although at some point it becomes impractical for the planar fit method 
where a matrix must be inverted to perform the multiple regression. 

The Lee (1998) approach uses regression of the vertical velocity on the 
horizontal wind speed for each wind direction category, generating a sequence 
of wind direction dependent coefficients a(<p) and The corrected mean 
vertical velocity (w) is then given by 

w = w m - ( a{<t > ) + &{<(>) U ) (2) 

where we use 36 wind direction categories of width 10 deg each and the mean 
horizontal wind speed is given by 

U = {u 2 + v 2 ) 1 /2 . (3) 

The planar fit approach uses multiple regression of the vertical velocity on 
the two components of the horizontal wind, and the corrected vertical motion 
is given by 


w = w m — (a + b u + c v) (4) 

where the coeffcients a, b and c are calculated using all the data. The tilt 
angle approach calculates a tilt angle as 
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T a = arctan(^p) (5) 

where W mn is the measured vertical velocity that has the experiment wide 
mean (potential instrument offset) removed. The procedure temporarily re- 
moves periods with the weakest mean flow ( U < 0.5 ms -1 ), where calculation 
of the tilt angle is not well posed and the tilt angle becomes erratic (e.g., Paw 
U et al., 2000). The corrected mean vertical velocity is then estimated for 
all the data, including the weakest wind periods, as 

w = w mn - U tan (T a (<£)). (6) 

where T a (<f>) is the bin-averaged wind direction dependent tilt angle using 36 
wind direction categories of width 10 deg. 

3.4 Results 

The three different tilt correction methods yield different w estimates (Figure 
2). Although the correlations between 30-minute average w estimates are 
high, the root-mean-square differences of a few cm s 1 are significant, and 
the maximum differences are as large as 10 cm s -1 . These large differences 
are discouraging for prospects of estimating mean vertical motion using a tilt 
correction method. 

The magnitude of the tilt angles found here (Figure 3) are comparable 
to or smaller than what we have typically found for other sites on relatively 
flat terrain with carefully levelled sensors. A sinusoidal dependence is a 
reasonable approximation to T a (<j>), implying that the planar fit method could 
be appropriate for this site. Differences between the bin-averaged T a (<f>) 
values and a pure sinusodial dependence, as implied by the planar fit method, 
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may be due to an insufficient number of samples (random sampling error) 
for the less frequent wind directions. On the other hand, some sites may be 
characterized by directionally dependent slopes that can not be reasonably 
represented by a planar fit, in which case the bin-averages of T a (<p) may be the 
better representation. In addition, higher order flow distortion effects may 
contribute to more complex azimuthal dependencies of the tilt angle, such as 
a double peak. For example, Hogstrom and Smedman (2004) highlight the 
calibration problems and flow distortion characteristics of Gill Solent sonic 
anemometers. Such problems may be less severe for the unobstructed design 
of the CSAT3. 

One cause of significant differences between the Lee and the tilt angle 
methods is demonstrated in Figure 4. For Lee’s method, the regression of 
w m on horizontal wind speed gives negative slope (negative fj) despite the 
fact that for 75% of the data, where the wind speed is less than 2 ms -1 , the 
measured average vertical velocity is positive and increases with increasing 
wind speed. The tilt angle method indicates very weak long-term mean 
positive vertical motion (T a (<f>) = 0.06 deg). For this wind direction category, 
the two diffferent approaches for correcting the vertical velocity generally 
have the opposite sign. Because Lee’s regression method includes a constant 
term a, the correction ( w — w m ) can be of either sign depending on wind 
speed, while the tilt angle method correction, which is near zero for this case, 
is always the same sign for a given wind direction. 

Lee’s method appears to fail for this example due to the combination of 
a nonlinear dependence of w m on wind speed and a strongly skewed wind 
speed distribution due to the high frequency of occurence of weak winds. The 
magnitude of the slope /3 in Lee’s method decreases with increasing averaging 
time by removal of some of the stronger wind speed cases. Recall that the 


12 


122 


vector average wind speed decreases with increasing averaging time due to 
wind direction variability. The tilt angle approach is more robust for this 
example in that the mean tilt angle is less sensitive to the choice of averaging 
time than are the coefficients from the regression approach. 

4 Divergence 

The horizontal coordinate system was rotated for each 100-s period such 
that horizontal u gradients are proportional to u at tower A minus u at 
tower C, and v gradients by tower E minus tower D. Horizontal gradients 
were evaluated over a fixed distance of approximately 200 m (Figure 1). 
No significant wind speed bias problems were found by an intercomparison 
study of the Handar anemometers. The horizontal gradients of u and v 
were calculated at the 3 measurement levels (1, 5 and 12 m) using finite 
differencing. Vertical integration of the horizontal gradients from the surface 
to 12 m was performed using a piecewise linear fit to the four levels, where 
the gradients are zero at z = 0. 

For the shorter towers D and E, where the measurements are limited to 
1 and 5 m, we assume the horizontal gradient of the mean wind is constant 
with height between 5 and 12 m. The measured gradients on the taller 
towers support this assumption. The fact that the horizontal gradient of 
the mean wind does not significantly increase above 5 m suggests that the 
primary mechanism generating divergence at this site is the surface roughness 
change. 

One issue with the divergence method of estimating w is the length scale 
over which the horizontal gradients are calculated. Ideally, the gradients 
would be calculated over a short distance in the immediate vicinity of the 
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vertical velocity measurements, however, the limited resolution of the hor- 
izontal wind speed measurements requires a sufficient seperation distance 
such that wind speed differences can be resolved. In some special circum- 
stances, such as with propogating divergence, the combination of seperation 
distance and averaging time could strongly influence the calculated diver- 
gence. However, if the divergence is primarily stationary and the horizontal 
wind components vary approximately linearily with distance, then the seper- 
ation distance and averaging time become less important. We did not find a 
strong averaging time dependence in the vertically integrated divergence. 

The divergence calculations indicate that this is a site of mean sinking 
motion for all wind directions and all times of day and night. We would 
expect to see weaker divergence due to the surface roughness change for those 
wind directions with longer fetch over the clearing prior to reaching the tower 
network. Theoretically, the divergence is a maximum at the upwind edge of 
the clearing and decreases with increasing distance (e.g., Garratt, 1990). The 
observations do indeed show that the minimum divergence occurs with north- 
west flow, where the tower network is located in the south-east corner of the 
clearing. 

5 Comparisons 

In this section we compare the w estimates calculated from w m and the 
different tilt correction methods (the direct methods) with the w estimate 
based on the horizontal divergence. There is large scatter for individual 
30- minute average estimates and small correlation (Figure 5). A 30- minute 
mean vertical motion of 10 cm s -1 , which implies an increase in the mean 
wind speed across the network of about 1.7 ms -1 , is probably too large to 
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be realistic and is not supported by the divergence measurements. Such 
magnitude vertical motions are found for all the direct methods, but not for 
the divergence method. This might indicate that the divergence approach is 
more robust. More spatial information is incorporated into the divergence 
calculation. If a 10 cm s -1 w was due to an error in specification of the tilt 
angle, the error in the ^dependence of T a would need to be about 3 deg 
for a mean wind speed of 2 ms -1 . It appears unlikely than an error of this 
magnitude occurs. 

The correlation between w from the tilt angle and planar fit methods 
and w from the divergence method increases with increasing averaging time, 
however, this is not the case for Lee’s method (Figure 6). The RMS w 
differences for all the direct methods decrease with increasing averaging time. 
This implies that a portion of the scatter is random. The standard deviation 
of the direct method w estimates is typically a factor of 2 greater than the 
standard deviation of W based on the divergence. 

All of the direct methods result in large scatter and large short-term fluc- 
tuations in w compared to the divergence method. Despite the apparent 
noise in the direct estimates, short-term fluctuations in w from the direct 
and divergence methods can be highly correlated (Figure 7). In this exam- 
ple, both the direct and divergence estimates clearly show enhanced sinking 
motion associated with stronger wind speeds, and enhanced rising motion 
with weaker winds. The direct method W fluctuations are too large to be due 
to errors in the tilt correction, despite being strongly wind speed dependent. 
The stronger amplitude response of the direct method w compared to the 
divergence estimate for this example remains unexplained. 

The agreement between the direct and divergence methods is worse when 
the magnitude of w from the direct methods does not increase monotonically 
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with height between the surface and 12 m (not shown). Such periods may 
represent complex flow situations that are impossible to resolve without a 
much finer grid of measurements. This situation occurred about one-half 
the time independent of the tilt correction method. Removing such periods 
improved the correlation between the 30-minute mean w estimates from the 
tilt angle and divergence methods from 0.39 to 0.47. These periods also 
include cases where w at both the 5 and 12-m levels is very small, in which 
case small random errors possibly due to flow distortion probably dominate 
the estimates. 

For very shallow nocturnal boundary layers, the vertical velocity mea- 
surement at 12 m may be partially decoupled from the horizontal flow field 
beneath it. To test this, we compared the estimates of vertical motion at 
5 m instead of 12 m above the surface. No significant improvement in the 
relationships was found at 5 m compared to 12 m (not shown). Interpre- 
tation is complicated by the fact that theoretically the relationship would 
be expected to degrade at 5 m compared to 12 m because the magnitude 
of uJ normally decreases with height (smaller signal to noise ratio) and any 
problems associated with surface heterogenity increase near the surface. 

As a further sensitivity study, the direct method w estimates were recal- 
culated at 12 m using the same procedures described above but only using 
nocturnal data. Despite this attempt to “tune” the tilt correction for noctur- 
nal conditions, the relationship between the direct vertical motion estimates 
and the divergence based estimates did not significantly improve (not shown). 
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5.1 Composites 

In the three month composite diurnal cycle, the W estimates from the pla- 
nar fit and tilt angle methods have some similarities to each other, and some 
similarity to the divergence based estimate of w (Figure 8). The similar com- 
posite diurnal behaviour is encouraging, although agreement on individual 
short time scales (30-minutes to 1-h) is the critical test for the ultimate goal 
of estimating vertical advection for process orientated studies. Contrary to 
the other methods tested, the Lee approach leads to significant mean rising 
motion in the afternoon. 

The direct and divergence methods disagree on the sign of the composite 
w during most of the night. Direct methods generally show rising motion of 
about 1 cm s _1 in contrast to weak sinking motion of 0.5 cm s _1 based on 
the divergence (Figure 8). The nocturnal sign difference could be due to the 
removal of the long-term mean w by the tilt correction methods. As discussed 
above (section 3), the tilt correction methods fail at sites with non-zero long- 
term mean vertical motion. Assuming that the divergence estimate is correct, 
and that long-term mean sinking motion due to the surface roughnesss change 
is realistic for this site, then the vertical advection of carbon dioxide that 
would be calculated using the direct method nocturnal w would be of the 
wrong sign for most of the night regardless of the choice of tilt correction 
method. 

The wind direction, wind speed and time of day dependencies of the com- 
posite w residuals for the tilt angle method and the uncorrected w m estimates 
are shown in Figure 9. The residual is defined as w for the direct method 
minus W from the divergence. The uncorrected estimate has the experiment 
wide mean (offset) of 0.9 cm s -1 removed, and not removing the offset would 
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make the uncorrected residuals in Figure 9 even larger positive. For the most 
common conditions of a 255 deg wind at 1.9 ms -1 , the composite residual 
is near zero for the tilt angle method. Assuming that the estimate based 
on the divergence is close to the true mean vertical motion, and comparing 
with the uncorrected estimate of vertical motion, it appears that the tilt an- 
gle method used here over-corrects in stronger winds and under-corrects in 
weaker winds. This is also seen in the diurnal pattern of the residual which 
is negative during the day, when the mean wind speed is significantly higher 
than at night. 

The positive residual for south-southwest flow (Figure 9) may be due in 
part to an exponential rather than linear fetch dependence of the divergence. 
An exponential dependence would be most pronounced near the upwind edge 
of the clearing, which is where the network is located for south-southwest 
flow. If the divergence decreased exponentially with fetch, the mean w corre- 
sponding to the divergence would be found in the upwind part of the network, 
not in the center of the network. As such, our direct measurement at the 
central tower site may underestimate the mean sinking motion. A larger neg- 
ative w from the divergence than from the direct method leads to a positive 
residual in Figure 9. 

6 Conclusions 

A three month field experiment was performed to contrast the mean vertical 
motion measured at a point on a central tower with the mean vertical motion 
calculated from the horizontal divergence measured using a network of towers 
surrounding the central tower. Three different tilt correction methods were 
applied. The RMS difference in the 30-minute average w estimates due to the 
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choice of tilt correction method was 2 to 3 cm s~ l . The 3-month composite 
diurnal pattern of w was sensitive to the choice of tilt correction method. 

The hour to hour variations in w from the direct methods are probably 
too large to be representative of a large spatial area. In such cases, the 
direct w estimates could be responding to divergence on horizontal scales 
smaller or larger than the divergence network. The RMS differences between 
the 30-minute average direct estimates of W and the estimates based on the 
divergence are 4 to 5 cm s -1 , and as a result, little confidence can be placed 
in the hour-to-hour variations in the w estimates, or in subsequent estimates 
of vertical advection. Of the three direct methods tested at this site, the tilt 
angle method best matched the divergence. 

The direct methods and the divergence method agree on some features 
of the composite diurnal cycle of w, however, the divergence indicates long- 
term mean sinking motion for all times of day and night, apparently due to 
the surface roughness change. Because the direct methods all remove the 
long-term mean, they fail to reproduce this result. As a consequence, the 
direct methods give the wrong sign of w at night, and thus would give the 
wrong sign of the nocturnal vertical advection of carbon dioxide. 

Our conclusion is that more confidence can be attached to vertical ad- 
vection estimates when the mean vertical motion is calculated from the hor- 
izontal divergence. A densor network of vertical velocity and divergence 
measurements would enable more definitive conclusions. Unfortuneatly, the 
measurements required to calculate the divergence are more costly and not 
normally available. 
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Figure 3: Mean and standard error of the tilt angle vs wind direction (100-s 
average data at 12 m). 
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Figure 4: Comparison of the Lee and tilt angle methods for wind direction 
category 220 < <j> < 230 deg. (100-s average data at 12 m). 
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Figure 5: 30-minute mean vertical motion at 12 m for the a) tilt angle, b) Lee, 
and c) planar fit methods vs the estimate based on the horizontal divergence. 
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Averaging Time (h) 

Figure 6: Averaging time dependence of the correlation and root mean square 
error for contrasting the three direct methods with the divergence method of 
estimating w. 






Hour Of Day 256 

Figure 7: Time series of 10-minute average mean wind speed divided by 10 
(dotted), w from direct tilt angle method (dashed), and w from divergence 
method (solid). 
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Figure 8: Mean and standard error of the three month composite diurnal 
cycle of mean vertical motion from four different methods. 
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Figure 9: Dependence of the composite w residual for the tilt angle method 
(solid), and for the uncorrected w m (dash), on wind direction, wind speed 
and time of day. 
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